SUBROUTINE simulatemun()

  USE Commonvars
  USE Random
  USE Statistics
  USE Nelder_Mead

  IMPLICIT NONE
  INCLUDE 'mpif.h'
  
  REAL(8), ALLOCATABLE :: simout(:,:,:), simoutt(:,:,:,:), simouttt(:,:,:,:), betart(:), betanort(:), notmayor_value_t(:,:,:), notmayor_value(:,:,:), ELNqcons_base_vec_t(:,:,:), ELNqcons_base_vec(:,:,:)
  REAL(8), EXTERNAL    :: zeroin, wealth_equivalent
  INTEGER, PARAMETER   :: N_other_var = 8
  REAL(8)              :: msavingf, stealf, qconsf, mwagef,  mnewsavingf, newstealf, newpconsf, newqconsf, yy, fundsf, utilf, ts1, ts2, ts3, zzpuf, zzprf, uud(6), pauditt(1:2), &
                          emaxr, emaxnor, numt(Nmomtemp), dent(Nmomtemp), num(Nmomtemp), den(Nmomtemp), frac_type_num(Ntypes,Nnterms), frac_type_numt(Ntypes,Nnterms), &
                          frac_type_dent(Nnterms), frac_type_den(Nnterms), frac_el_app_num(Nrcont,Nnterms), frac_el_app_numt(Nrcont,Nnterms), frac_el_app_dent(Nnterms), frac_el_app_den(Nnterms), &
                          frac_type_el_app_num(Ntypes,Nrcont,Nnterms), frac_type_el_app_numt(Ntypes,Nrcont,Nnterms), frac_type_el_app_dent(Nnterms), frac_type_el_app_den(Nnterms), &
                          frac_type_run_num(Ntypes), frac_type_run_numt(Ntypes), frac_type_run_den(1), frac_type_run_dent(1), other_var_numt(N_other_var), other_var_dent(N_other_var), &
                          other_var_num(N_other_var), other_var_den(N_other_var), other_var(N_other_var), fr_stolen, dum_steal, dum_steal_x, vemaxpf, &
                          ln_relqcons, qconsvec(Nq), d_run, fine_fun, rundec_netW, v_slope, multi_slope, zzpr_shock_t, run_shock_t, thres_frac_stolen, newstealf_nome, &
                          ability, potential_wage, wage_pastmayor, msaving_bf, med_wealth(3), temp(2*Nsimterms), frac_fine_p, frac_fine_c, &
                          wealth_star, a, b, fa, fb, fr_wealth, will_to_pay, vemax_ini, vemax_base, vemax_policy, av_will_to_pay_t(Nactmun,Nsimterms), av_will_to_pay(Nactmun,Nsimterms), &
                          notmayor_base_value, ELNqcons_base, ELNqcons, temp_var(Nvar), var_base(Nsim,Nactmun,Nsimterms,2), outelect_base, msaving_base, &
                          ability_unit, wealth_unit, Pr_i_win_temp, Pwin(2)
  INTEGER              :: i, j, index_steal_x, tt, outelect, auditf,auditf3, thief3, mayorage, mayorage_elec, incumb_age, incumb_educ, nterms, ipastf, Totterms, &
                          aenoc_mun, aud, aud3, ipau, ipauf, totloop, ntasks, restasks, counter, iw, is, ifu, ipr, ist, imun, ieducf, isav, iab, infeasible, dft, &
                          nsimt, thief, ub_sav, ub_sav_bi, ub_ab_bi, el_ap, med_edu(3), med_age(3), med_ipast(3), med_edu_i, med_ipast_i

  EXTERNAL             :: objfcn_ml, notmayordecision, rmultinomial, mayordecision

  INTERFACE
     SUBROUTINE sfun(n, x, f, g)
       IMPLICIT NONE
       INTEGER, PARAMETER      :: dp = SELECTED_REAL_KIND(12, 60)
       INTEGER, INTENT(IN)     :: n
       REAL (dp), INTENT(IN)   :: x(:)
       REAL (dp), INTENT(OUT)  :: f, g(:)
     END SUBROUTINE sfun
     SUBROUTINE sfun1(n, x, f, g)
       IMPLICIT NONE
       INTEGER, PARAMETER      :: dp = SELECTED_REAL_KIND(12, 60)
       INTEGER, INTENT(IN)     :: n
       REAL (dp), INTENT(IN)   :: x(:)
       REAL (dp), INTENT(OUT)  :: f, g(:)
     END SUBROUTINE sfun1
  END INTERFACE

  if (myrank == 0) write(*,*) 'here 2'
  
  CALL cpu_time(ts1)

  flag_simulation = 1

  !Compute the size of all loop
  totloop = Nsim
  ntasks = INT(totloop/nprocs)
  restasks = totloop - ntasks*nprocs

  ALLOCATE(simout(Nactmun,Nsimterms,Nvar))
  simout = 0d0

  IF (flag_simulate_data == 1) THEN
     ALLOCATE(simoutt(Nactmun,Nsim,Nsimterms,Nvar))
     simoutt = 0d0
  END IF

  ALLOCATE(notmayor_value_t(Nactmun,Nsim,Nsimterms), notmayor_value(Nactmun,Nsim,Nsimterms), ELNqcons_base_vec_t(Nactmun,Nsim,Nsimterms), ELNqcons_base_vec(Nactmun,Nsim,Nsimterms))
  notmayor_value_t = 0d0
  notmayor_value = 0d0
  ELNqcons_base_vec_t = 0d0
  ELNqcons_base_vec = 0d0

  av_will_to_pay_t = 0d0
  !In 2019 in Brazil, median age is 33.2 (2=32-37 in our code) and average education is 9.3 (1=high school), We use median age, education, and wealth conditional on the size of the municipality, computed in structest_Feb_2020.do
  !We allow the median value to vary with municipality size
  med_edu(1:2) = 1
  med_edu(3) = 1
  med_age(1) = 2
  med_age(2:3) = 2
  med_ipast(:) = 1
!!$  Median and mean wealth in Brazil:
!!$  The numbers are in US dollars, we have to multiply by X to get number in Brazilian Reals; What exchange rate do we use? The current one is $1 = 5.2335 Reals, the average for 2019 is $1 = 3.9457 Reals
!!$  for 2001-2010 (mean) for 2010 (median): https://www.credit-suisse.com/media/assets/corporate/docs/about-us/research/publications/credit-suisse-global-wealth-databook.pdf
!!$  For 2019: https://www.credit-suisse.com/media/assets/corporate/docs/about-us/research/publications/global-wealth-databook-2019.pdf
!!$  for 2019: https://en.wikipedia.org/wiki/List_of_countries_by_wealth_per_adult
!!$  med_wealth(1) = 218.248d0
!!$  med_wealth(2) = 287.3192d0
!!$  med_wealth(3) = 430.7853d0
  med_wealth(1) = 47.38653233d0 ! = 14,242 * 0.633892551 * 5.2489 / norm (average per-adult brazilian wealth in US dollars in 2006 in Table 2.4 times price index to deflate to 2000 times the average exchange rate in April 1st 2020) !66.9727943d0 ! = 24,289 * 0.525315827 * 5.2489 / norm (average per-adult brazilian wealth in US dollars in 2010 in Table 2.4 times price index to deflate to 2000 times the average exchange rate in April 1st 2020) ! 92.921235d0 ! = 23,550 * 3.9457/norm (average per-adult brazilian wealth in US dollars in 2019 in Table 21 times the average exchange rate in 2019).
  med_wealth(2) = 47.38653233d0 !66.9727943d0 ! = 24,289 * 0.525315827 * 5.2489 / norm (average per-adult brazilian wealth in US dollars in 2010 in Table 2.4 times price index to deflate to 2000 times the average exchange rate in April 1st 2020).
  med_wealth(3) = 47.38653233d0 !66.9727943d0
  !Given the choosen wealth, education, and age, total initial resources in the computation of wtp = 156.272929 = 92.921235d0*1.05+pmwage*12
  !Wtp as fraction of wealth should therefore be computed as wtp/156.272929
  !We read in the residents' values for the base case
  IF (flag_simulate_data == 1 .AND. flag_base == 0) THEN
     OPEN(unit=28, FILE="/home/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/Values_Base.dat")
     DO j = 1, Nsim
        DO i = 1, Nactmun
           READ(28,*) temp(:)
           notmayor_value(i,j,:) = temp(1:Nsimterms)
           ELNqcons_base_vec(i,j,:) = temp(Nsimterms+1:2*Nsimterms)
!!$           IF (myrank == 0) write(*,'(a,2i6,6f20.10)') 'notmayor_value(i,j,:), ELNqcons_base_vec(i,j,:)', j, i, notmayor_value(i,j,:), ELNqcons_base_vec(i,j,:)
        END DO
     END DO
     CLOSE(28)
  END IF

  Totterms = Nsimterms
  numt = 0d0
  dent = 0d0
  num = 0d0
  den = 0d0
  
  frac_type_numt = 0d0
  frac_type_dent = 0d0
  frac_type_num = 0d0
  frac_type_den = 0d0

  frac_type_run_numt = 0d0
  frac_type_run_dent = 0d0
  frac_type_run_num = 0d0
  frac_type_run_den = 0d0

  frac_el_app_numt = 0d0
  frac_el_app_dent = 0d0
  frac_el_app_num = 0d0
  frac_el_app_den = 0d0
  
  frac_type_el_app_numt = 0d0
  frac_type_el_app_dent = 0d0
  frac_type_el_app_num = 0d0
  frac_type_el_app_den = 0d0

  other_var_numt = 0d0
  other_var_dent = 0d0
  other_var_num = 0d0
  other_var_den = 0d0

  frac_fine_c = 0d0

!!$        Nsimt = 1
  Nsimt = nsim

  ! The following integers are only used in the computation of emax and not in the simulation,
  ! but we need them in mayordecision.f90
  iw = 0
  is = 0
  ifu = 0
  ipr = 0

  ! If we control for observable heterogeneity at the mayor level
  IF (flag_no_selec_obs == 1) THEN

     OPEN(unit=18, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_base_Jan2021.dat")  

     DO j = 1, Nsim
        DO i = 1, Nactmun
           DO tt = 1, Nsimterms
              
              READ(18,*) temp_var
              ! outelect base
              var_base(j,i,tt,1) = temp_var(15)
              ! savings base, we have to use var 12 (savings entering the period) and not var 16 (optimal saving at the end)
              var_base(j,i,tt,2) = temp_var(12)
              
           END DO
        END DO
     ENDDO
     
  END IF

  counter = 0
  DO j = 1, Nsimt

     counter = counter + 1
     IF ((counter <= ntasks*nprocs .AND. counter >= myrank*ntasks + 1 .AND. counter <= myrank*ntasks + ntasks) .OR. &
          (counter > ntasks*nprocs .AND. counter - ntasks*nprocs - 1 == myrank)) THEN
        
        DO i = 1, Nactmun

           CALL cpu_time(ts3)
           
           DO tt = 1, Totterms

              infeasible = 0
              
              thief = 0
              thief3 = 0
              aud3 = 0 
              auditf3 = 0

              !We multiply the zzpr_shocks distributed N(0,1) times the std deviation of the distribution to get a N(0,sigma2)
              zzpr_shock_t = sigma_zzpr*zzpr_shock(i,j,tt)
!!$              zzpr_shock_t = 0d0
              !Shock to the decision to run
              run_shock_t = sigma_run*rnormprob(i,j,tt) !zzpr_shock(i,j,tt)
!!$              if (myrank == 0) write(*,*) 'run_shock_t', run_shock_t, sigma_run, rnormprob(i,j,tt)

              !We transform the draw from the standard normal in draws from a log-normal(mu_fine,sigma_fine)
              !Fraction of the fine applied in case the mayor is audited this term
              frac_fine_p = frac_fine_c
              !Fraction of the fine used in mayordecisions to determine the amount of wealth he will have in the next term if he steals
              !frac_fine_p must always be equal frac_fine_n in the previous term: what the mayor believes will be applied if audited, must be equal to what will actually be applied if audited 
              frac_fine_c = EXP(sigma_fine*fine_norm(i,j,tt) + mu_fine) 
!!$              frac_fine_c = EXP(sigma_fine*fine_norm(i,1,tt) + mu_fine) 
              IF (frac_fine_c > max_frac_fine) frac_fine_c = max_frac_fine

              uud(1:6) = udraw(i,j,tt,:)

              IF (flag_no_selec_obs == 1) THEN
                 ! electoral outcome in the base case at the start of period tt
                 outelect_base = var_base(j,i,tt,1)
                 ! savings in the base case at the beginning of the period
                 msaving_base = var_base(j,i,tt,2)
              END IF

              IF (tt == 1) THEN

                 !In the first term set the prob of being audited to 0.04 following the data in Brazil
                 ipau = 1
                 ipauf = 1
                 ! In tt = 1, all mayors are term 1, so have no past
                 ipastf = 1
                 ! Since they have no past, their past stealing is equalt to zero
                 stealf = 0d0
                 ! And the previous incumbents decided not to run
                 d_run = 0d0

                 !Observed/Conditioning Variables
                 !imun descibes the type of municipality in terms of population; we have six types.
                 imun = INT(initdata(i,10))
                 popsizeg = popsize(imun)
                 auditf = INT(initdata(i,2))
                 mayorage = INT(initdata(i,3))
                 ! For now, censoring the age of mayors: later we clean up the sample
                 mayorage_elec = MIN(Retage-1,mayorage)
                 incumb_age = mayorage_elec
                 ieducf = INT(initdata(i,4))

                 fundsf = per_funds*initdata(i,6)/norm
                 ! Since they have no past, their past fraction stolen is equalt to zero
                 fr_stolen = stealf/fundsf

!!$                 zzprf = initdata(i,5) * EXP(zzpr_shock_t)
                 zzprf = initdata(i,5)
                 nterms = INT(initdata(i,8))

                 !We don't need qconsf for Brasil since all mayors are in their first term in 1996
                 !However we need to give it a value for drawinitunob

                 IF (imun == 1) THEN
                    qconsvec = qconsvec1
                 ELSE IF (imun == 2) THEN
                    qconsvec = qconsvec2
                 ELSE IF (imun == 3) THEN
                    qconsvec = qconsvec3
                 END IF
                 !We compute EXP(qconsvec) to have the level
                 qconsf = EXP(qconsvec(2))

                 !Wage of the mayor, it depends on the type of municipality
                 !We use the wage from our grid                 
!!$                 mwagef =  n_times_wage*initdata(i,7)/norm
                 mwagef =  wage_mayor(imun)

                 ! Potential wage in the private sector for the incumbent
                 potential_wage = wage_pastmayor(mayorage,imun,ieducf,ipastf)

                 !Outcome of the election in the first period, 1 means that the incumbent wins the election

                 !In the data, at tt = 1, all mayors are new mayors
                 outelect = INT(initdata(i,9))       
                 yy = 1d0

                 IF (outelect == 0 .AND. auditf == 1) WRITE(*,*) 'wrong initial data, outelect, auditf',  outelect, auditf

                 !We draw ability from a log-normal distribution
!!$                 ability = EXP(sigma_a*ability_norm(i,1,tt)+mu_a)
                 ability = EXP(sigma_a*ability_norm(i,j,tt)+mu_a)
!!$                 ability = sigma_a*ability_norm(i,j,tt)+mu_a

                 !If we fix ability at some level
                 IF (flag_fixed_ability == 1) ability = fixed_ability

!!$                 IF (myrank == 0) write(*,'(a,3i5,4f20.10)') 'ability inc 1',i,j,tt, ability, sigma_a, ability_norm(i,j,tt), mu_a

                 !We draw the electoral appeal type
                 ! 1 is low electoral appeal, 2 is high electoral appeal, pappeal(1) is the probability of a low electoral appeal type
                 el_ap = 1
                 IF (uud(6) >= pappeal(1)) THEN
                    el_ap = 2
                 END IF

!!$                 if (myrank == 0) write(*,'(a,1i4,2f20.10)') 'el_ap', el_ap, uud(6), pappeal(1)

                 msavingf = initdata(i,18)/norm

!!$                 ! If we simulated the base case without selection on observables in the first and second terms, then flag_no_selec_obs = 1 and flag_no_selec_obs_from_1st = 1, and we assume that all mayors have the average variables observed in the data that are not mayor specific; if flag_no_selec_obs_from_1st == 1 then we set all observables equal to the average or median
!!$                 IF (flag_no_selec_obs == 1 .AND. flag_no_selec_mun_obs_from_1st == 1) THEN
!!$                    IF (flag_no_selec_obs_from_1st == 1) THEN
!!$                       mayorage = med_age(imun)
!!$                       mayorage_elec = MIN(Retage-1,mayorage)
!!$                       ieducf = med_edu(imun)
!!$                       msavingf = av_saving
!!$                    END IF
!!$                    zzprf = av_zzpr
!!$                    fundsf = av_funds
!!$                 END IF

!!$                 if (myrank == 0) write(*,'(a,2f20.10)') 'savings', msavingf, initdata(i,18)

!!$                 !We draw savings from a non-parameteric distribution because we don't know them in the main sample
!!$                 IF (imun == 1) THEN
!!$                    Nmunsav = Ntotmunnoc_mun1
!!$                    ALLOCATE(nonpardist_sav(Nmunsav,3),nonpardist_sav_temp(Nmunsav))
!!$                    !age
!!$                    nonpardist_sav(:,1) = nonpardistnoc_mun1(:,1,1)
!!$                    !education
!!$                    nonpardist_sav(:,2) = nonpardistnoc_mun1(:,1,2)
!!$                    !wealth
!!$                    nonpardist_sav(:,3) = nonpardistnoc_mun1(:,1,5)
!!$                    nonpardist_sav_temp = 0d0
!!$                 ELSE IF (imun == 2) THEN
!!$                    Nmunsav = Ntotmunnoc_mun2
!!$                    ALLOCATE(nonpardist_sav(Nmunsav,3),nonpardist_sav_temp(Nmunsav))
!!$                    !age
!!$                    nonpardist_sav(:,1) = nonpardistnoc_mun2(:,1,1)
!!$                    !education
!!$                    nonpardist_sav(:,2) = nonpardistnoc_mun2(:,1,2)
!!$                    !wealth
!!$                    nonpardist_sav(:,3) = nonpardistnoc_mun2(:,1,5)
!!$                    nonpardist_sav_temp = 0d0
!!$                 ELSE IF (imun == 3) THEN
!!$                    Nmunsav = Ntotmunnoc_mun3
!!$                    ALLOCATE(nonpardist_sav(Nmunsav,3),nonpardist_sav_temp(Nmunsav))
!!$                    !age
!!$                    nonpardist_sav(:,1) = nonpardistnoc_mun3(:,1,1)
!!$                    !education
!!$                    nonpardist_sav(:,2) = nonpardistnoc_mun3(:,1,2)
!!$                    !wealth
!!$                    nonpardist_sav(:,3) = nonpardistnoc_mun3(:,1,5)
!!$                    nonpardist_sav_temp = 0d0
!!$                 END IF
!!$                 
!!$                 ii = 1
!!$                 DO jj = 1, Nmunsav
!!$                    IF (INT(nonpardist_sav(jj,1)) == mayorage_elec .AND. INT(nonpardist_sav(jj,2)) == ieducf) THEN
!!$                       
!!$                       nonpardist_sav_temp(ii) = nonpardist_sav(jj,3)
!!$                       ii = ii + 1
!!$                       
!!$                    END IF
!!$                 END DO
!!$
!!$                 aenoc = INT(uud(3) * (DBLE(ii)-1.0d0)) + 1
!!$                 msavingf = nonpardist_sav_temp(aenoc)
!!$
!!$                 DEALLOCATE(nonpardist_sav,nonpardist_sav_temp)
                 
!!$                 write(*,*) 'aenoc', myrank, aenoc, msavingf

                 !In the first period this is irrelevant because all mayors are in their 1st term
                 thief=0
!!$                 IF(stealf > l_steal_prob/norm) thief =1
                 IF(stealf/fundsf > l_fr_stolen) thief =1

                 !We save the variables for the current period
                 simout(i,tt,1) = DBLE(i) 
                 simout(i,tt,2) = DBLE(j) 
                 simout(i,tt,3) = DBLE(tt) 
                 simout(i,tt,4) = fr_stolen
                 !qconsf is in levels
                 simout(i,tt,5) = qconsf
                 simout(i,tt,6) = mwagef
                 simout(i,tt,7) = DBLE(auditf)
                 simout(i,tt,8) = elap_vec(el_ap)
                 simout(i,tt,9) = fundsf
!!$                 simout(i,tt,10) = zzprf/EXP(zzpr_shock_t)
                 simout(i,tt,10) = zzprf
                 simout(i,tt,11) = ability
                 
                 !This is only useful is we write out data
                 IF (flag_simulate_data == 1) THEN
                    simoutt(i,j,tt,1) = DBLE(i) 
                    simoutt(i,j,tt,2) = DBLE(j) 
                    simoutt(i,j,tt,3) = DBLE(tt) 
                    simoutt(i,j,tt,4) = fr_stolen
                    simoutt(i,j,tt,5) = qconsf
                    simoutt(i,j,tt,6) = mwagef
                    simoutt(i,j,tt,7) = DBLE(auditf)
                    simoutt(i,j,tt,8) = elap_vec(el_ap)
                    simoutt(i,j,tt,9) = fundsf
!!$                    simoutt(i,j,tt,10) = zzprf/EXP(zzpr_shock_t)
                    simoutt(i,j,tt,10) = zzprf
                    simoutt(i,j,tt,11) = ability
                 END IF

              ELSE !IF tt==1
                 
                 !In the first term set the prob of being audited to 0.04 following the data in Brazil
                 IF (tt == 2) THEN
                    ipau = 1
                    ipauf = 2
                 ELSE
                    ipau = 2
                    ipauf = 2                    
                 END IF
                 
                 !If nterms >= Nnterms the mayor must step down
                 dft = 0
                 IF (nterms >= Nnterms) dft = 1
 
                 !Incumbent starting variables
                 stealf = newstealf
                 ! Both stealf and funds are for the previous period
                 fr_stolen = stealf/fundsf
                 !The following index is needed to choose the betas for a mayor who stole more than x%
                 index_steal_x = 1
                 IF (fr_stolen > l_fr_stolen_x) index_steal_x = 2                 
                 !qconsf is in levels
                 qconsf = newqconsf
                 msavingf = mnewsavingf                 
!!$                 IF (flag_no_selec_obs == 1 .AND. flag_no_selec_sav == 1) msavingf = msaving_base

                 ieducf = ieducf
                 incumb_educ = ieducf
                 mayorage = mayorage_elec + 1
                 incumb_age = mayorage

!!$                 IF (flag_no_selec_sav == 1) msavingf = av_saving
!!$                 ! If we simulated the base case without selection on observables in the second term, then flag_no_selec_obs = 1 and we do not update the mayor's observables, but we update the choices (variables above)
!!$                 IF (flag_no_selec_obs_from_1st == 0) THEN
!!$                    ieducf = ieducf
!!$                    mayorage = mayorage_elec + 1
!!$                    incumb_age = mayorage
!!$                 END IF

!!$                 write(*,'(a,5i4,2f20.10)') 'outelect_base 0', j, i, tt, mayorage, ieducf, outelect_base, msavingf

                 !audit
                 pauditt(1:2) = paudit(0:1,ipau)
                 IF (flag_audit_nterm == 1 .AND. (nterms == Nnterms .OR. mayorage >= Retage)) pauditt(1:2) = paudit_nterm(0:1)
                 CALL rmultinomial(uud(5), 2, pauditt, aud)
                 auditf = aud-1

                 !Policy that audits all mayors that were audited and caught stealing in the past
                 !This is the right place for setting auditf = 1, beacause the policy has an effect only if the mayor was caught stealing in the past and not in the current period
                 !Notice that the policy has an effect only through the value function if Nnterms<=2, because only mayors at the end of the second term can have ipastf = 3
                 !If Nnterms>2, then it can also have a direct effect
                 IF (flag_policy_audit == 1) THEN
                    IF (ipastf == 3) auditf = 1
                 END IF

                 !If the mayor has served the maximum number of terms or s/he is older than the maximum age (Retage)
                 !the challenger automatically wins the election.
                 IF (dft == 1 .OR. mayorage >= Retage .OR. infeasible == 1) THEN
                    yy = -99999999999d0
                    outelect = -999999
                    infeasible = 0
                    d_run = -1d0
                    GO TO 80
                 END IF

                 IF (infeasible == 1) THEN
                    yy = -77777777777777d0
                    outelect = -777777
                    infeasible = 0
                    d_run = -777d0
                    GO TO 80
                 END IF

                 !ln_relqcons is first used in the computation of the value function
                 !the value function depends on log of public consumption acounting for the degree of rivalry: qconsvec = log(Q/P**eta)
                 ln_relqcons = LOG(qconsf)-eta*LOG(popsizeg)

                 !The first decision of the incumbent is whether to run or not.
                 !We first compute the mayor's value function if s/he decides to run this period 
                 !ALL THIS IS CONDITIONAL ON KNOWING THE RESULT OF AUDIT!
                 !For now, this is only valid with Nnterms == 2; it has to be modified to allow for iterm =3 or more!

!!$                 IF (msavingf < 0d0) write(*,'(a,1f20.10)') 'neg sav 1 1', msavingf

                 msaving_bf = msavingf
                 IF(auditf == 0) THEN
                    !In this case, we keep the ipastf from the previous term to use the correct wage for the tax policy wage
                    rundec_netW = msavingf
!!$                 ELSE IF(auditf == 1 .AND. stealf <= l_steal_prob/norm) THEN
                 ! fr_stolen is for the previous period
                 ELSE IF(auditf == 1 .AND. fr_stolen <= l_fr_stolen) THEN
                    ipastf = 1
                    IF (Npast > 1) THEN
                       !If ipastf == 3, we keep the value from the previous term to use the correct wage for the tax policy wage
                       IF (ipastf < 3) THEN
                          ipastf = 2
                       END IF
                    END IF
                    rundec_netW = msavingf
!!$                 ELSE IF(auditf == 1 .AND. stealf > l_steal_prob/norm) THEN
                 ! fr_stolen is for the previous period
                 ELSE IF(auditf == 1 .AND. fr_stolen > l_fr_stolen) THEN
                    !We have to use the frac_fine saved from the previous term
                    frac_fine_g = frac_fine_p
                    rundec_netW = msavingf - fine_fun(stealf)
                    msavingf = rundec_netW
                    ipastf = 1
                    IF (Npast > 1) THEN
                       ipastf = 3
                    END IF
                 END IF

!!$                 ! We add this if we want previous stealing to have no effect on savings when we set savings = average savings
!!$                 IF (flag_no_selec_sav == 1) msavingf = av_saving

                 ! Potential wage in the private sector for the incumbent
                 potential_wage = wage_pastmayor(mayorage,imun,ieducf,ipastf)

                 ALLOCATE(betart(MMnew))
                 betart =  betasimrun(mayorage,imun,el_ap,ipastf,nterms+1,ipauf,ieducf,:)

                 !Find the point on the saving domain where we need to compute the approximated value function
                 DO isav = 1, Nwealth
                    IF (rundec_netW <  wealth(isav)) THEN
                       ub_sav = isav
                       EXIT
                    END IF
                    ub_sav = Nwealth
                 END DO
                       
                 !COMPUTE THE LINEAR APPROXIMATION
                 IF (ub_sav == 1) THEN
                          
                    multi_slope = DBLE(CEILING(ABS((rundec_netW - wealth(1)))/inc_wealth))
                    V_slope = (1d0+multi_slope)*(betart(2) - betart(1))/(wealth(2) - wealth(1))
                    
                 ELSE
                    
                    V_slope = (betart(ub_sav) - betart(ub_sav-1))/(wealth(ub_sav) - wealth(ub_sav-1))
                    
                 END IF
                 
                 emaxr = betart(MMnew-Npol_ab+1) + betart(MMnew-Npol_ab+2)*ability + betart(MMnew-Npol_ab+3)*ability**2d0 + V_slope*(rundec_netW - wealth(ub_sav)) + betart(ub_sav)

!!$                 if (myrank == 0) write(*,*) 'emaxr',i,j,tt, emaxr, betart(MMnew-Npol_ab+1), betart(MMnew-Npol_ab+2)*ability, betart(MMnew-Npol_ab+3)*ability**2d0, V_slope*(rundec_netW - wealth(ub_sav)), betart(ub_sav), run_shock_t
!!$                 if (myrank == 0) write(*,*) 'betart', i,j,tt,mayorage,nterms,betart
!!$                 if (myrank == 0) write(*,'(a,1i5,7f20.10)') 'emaxr', ub_sav, V_slope, (rundec_netW - wealth(ub_sav)), rundec_netW, wealth(ub_sav), betart(ub_sav), betart(ub_sav-1)

                 !We add a shock to the value to run
                 emaxr = emaxr + run_shock_t
                 ! If we evaluate the no-run policy it is easier to analyse the data without the shock
                 IF (flag_no_run == 1 .AND. auditf == 1 .AND. fr_stolen > l_fr_stolen) emaxr = emaxr - run_shock_t
                 IF (flag_no_run_x == 1 .AND. auditf == 1 .AND. fr_stolen > l_fr_stolen_x) emaxr = emaxr - run_shock_t

                 !We then compute the mayor's value function if s/he decides not to run this period 
                 !This is the expected value for this period. No exogenous variables (other than audit) is realized
                 ALLOCATE(betanort(Nwealth))
                 betanort = betasimnor(mayorage,imun,ipastf,ipauf,ieducf,:)

                 !COMPUTE THE LINEAR APPROXIMATION
                 IF (ub_sav == 1) THEN
                    
                    multi_slope = DBLE(CEILING(ABS((rundec_netW - wealth(1)))/inc_wealth))
                    V_slope = (1d0+multi_slope)*(betanort(2) - betanort(1))/(wealth(2) - wealth(1))
                    
                 ELSE
                    
                    V_slope = (betanort(ub_sav) - betanort(ub_sav-1))/(wealth(ub_sav) - wealth(ub_sav-1))
                    
                 END IF
                 
                 emaxnor = V_slope*(rundec_netW - wealth(ub_sav)) + betanort(ub_sav)

!!$                 if (fr_stolen>0d0 .and. auditf>0) THEN
!!$                    if (myrank == 0) write(*,'(a,9f20.10,2i5)') 'emaxr nor', emaxr-run_shock_t, emaxnor, emaxr, run_shock_t, fr_stolen, ability, ln_relqcons, rundec_netW, msaving_bf, ub_sav, auditf
!!$                 end if

!!$                 DEALLOCATE(betart,betanort)

                 IF (emaxnor > emaxr) THEN
                    yy = -111111111d0
                    outelect = -111111
                    !Dummy to record whether the incumbent run for reelection
                    d_run = 0d0

!!$                    if (myrank == 0 .and. tt == 3 .and. auditf == 1) THEN
!!$                       write(*,'(a,1i5,5f20.10,1i5)') 'NOT RUN', tt, emaxr-run_shock_t, emaxnor, emaxr, run_shock_t, fr_stolen, ub_sav
!!$                    end if

                 ELSE
                    
                    !Dummy to record whether the incumbent run for reelection
                    d_run = 1d0

                    ! Find the point on the ability grid that corresponds to the actual abilty; we need it because the probability of winning the election is defined for discrete values of wealth and ability
                    ! Ability is log-normally distributed, see momentsc.f90 line 244
                    ! ub_ab is used inside Pr_i_win (in Pr_i_win ub_ab_app cannot be above Nab)
                    DO iab = 2, Nab
                       IF (ability <  a_discr(iab)) THEN
                          ub_ab_bi = iab
                          EXIT
                       END IF
                       ub_ab_bi = Nab
                    END DO
                    
                    ! ind the point on the wealth grid that corresponds to the actual wealth;
                    DO isav = 2, Nwealth
                       IF (rundec_netW <  wealth(isav)) THEN
                          ub_sav_bi = isav
                          EXIT
                       END IF
                       ub_sav_bi = Nwealth
                    END DO
                 
                    ! We only know the reelection probability for discrete value of ability and wealth
                    ! Given the actual value of ability and wealth, we know in which interval of ability and wealth they reside given our grids
                    ! We therefore know the rectangle in which the point (ability,wealth) resides: [ub_sav-1, ub_sav; ub_ab-1, ub_ab]
                    ! We also know the reelection probability at each point of the rectangle
                    ! We can then use bilinear interpolation to approximate the reelection probability at the point (ability,wealth) using the probabilities at each corner
                    ! See https://blogs.sas.com/content/iml/2020/05/18/what-is-bilinear-interpolation.html for the formulas
                    ! We use ability on the horizontal axis and wealth on the vertical axes, hence
                    ! lower-left corner = ub_ab-1, ub_sav-1, upper-left corner = ub_ab-1, ub_sav, lower-right corner = ub_ab, ub_sav-1, upper-right corner = ub_ab, ub_sav
                    
                    ! 1st step: map our rectangle onto the unit square:
                    IF (ability >= a_discr(1) .AND. ability <= a_discr(Nab)) THEN
                       ability_unit = (ability - a_discr(ub_ab_bi-1))/(a_discr(ub_ab_bi) - a_discr(ub_ab_bi-1))
                    ELSE IF (ability < a_discr(1)) THEN
                       ability_unit = 0d0
                    ELSE IF (ability > a_discr(Nab)) THEN
                       ability_unit = 1d0
                    END IF
                    
                    IF (rundec_netW >= wealth(1) .AND. rundec_netW <= wealth(Nwealth)) THEN
                       wealth_unit = (rundec_netW - wealth(ub_sav_bi-1))/(wealth(ub_sav_bi) - wealth(ub_sav_bi-1))
                    ELSE IF (rundec_netW < wealth(1)) THEN
                       wealth_unit = 0d0
                    ELSE IF (rundec_netW > wealth(Nwealth)) THEN
                       wealth_unit = 1d0
                    END IF
                    
                    ! 2nd step: we do the bilinear interpolation in the unit square:
                    Pr_i_win_temp = Pr_i_win(imun,mayorage,ipastf,el_ap,ieducf,ub_sav_bi-1,ub_ab_bi-1,ipauf,nterms)*(1d0-ability_unit)*(1d0-wealth_unit) + &
                         Pr_i_win(imun,mayorage,ipastf,el_ap,ieducf,ub_sav_bi-1,ub_ab_bi,ipauf,nterms)*ability_unit*(1d0-wealth_unit) + &
                         Pr_i_win(imun,mayorage,ipastf,el_ap,ieducf,ub_sav_bi,ub_ab_bi-1,ipauf,nterms)*(1d0-ability_unit)*wealth_unit + &
                         Pr_i_win(imun,mayorage,ipastf,el_ap,ieducf,ub_sav_bi,ub_ab_bi,ipauf,nterms)*ability_unit*wealth_unit

                    Pwin = (/ 1d0 - Pr_i_win_temp, Pr_i_win_temp /)
                    CALL rmultinomial(uud(3), 2, Pwin, outelect)
                    outelect = outelect - 1
                    yy = 1d0
                    IF (outelect == 0) yy = 01d0

                    IF (Pr_i_win_temp <= 0d0) THEN

                       IF (myrank == 0) write(*,'(a,2i4,5f20.10)') 'unit variables', ub_ab_bi, ub_sav_bi, ability_unit, wealth_unit, rundec_netW, wealth(ub_sav_bi-1), wealth(ub_sav_bi)
                       
                       IF (myrank == 0) write(*,'(a,9i4)') 'indeces', imun,mayorage,ipastf,el_ap,ieducf,ub_sav_bi-1,ub_ab_bi-1,ipauf,nterms
                       
                       IF (myrank == 0) write(*,'(a,9f20.10)') 'Pr_i_win_temp', Pr_i_win_temp, Pr_i_win(imun,mayorage,ipastf,el_ap,ieducf,ub_sav_bi-1,ub_ab_bi-1,ipauf,nterms), (1d0-wealth_unit), &
                            Pr_i_win(imun,mayorage,ipastf,el_ap,ieducf,ub_sav_bi-1,ub_ab_bi,ipauf,nterms), ability_unit*(1d0-wealth_unit), &
                            Pr_i_win(imun,mayorage,ipastf,el_ap,ieducf,ub_sav_bi,ub_ab_bi-1,ipauf,nterms), (1d0-ability_unit)*wealth_unit, &
                            Pr_i_win(imun,mayorage,ipastf,el_ap,ieducf,ub_sav_bi,ub_ab_bi,ipauf,nterms), ability_unit*wealth_unit
                       
                       IF (myrank == 0) write(*,'(a,1i4,2f20.10)') 'outelect', outelect, uud(3), 1d0 - Pr_i_win_temp

                    END IF
                    
                    dum_steal = 0d0
                    dum_steal_x = 0d0
                    IF (fr_stolen > l_fr_stolen) dum_steal = 1d0
                    IF (fr_stolen > l_fr_stolen_x) dum_steal_x = 1d0
                    
                    ! If we want to implement a policy in which the mayor looses the election with probability 1 if he steals (equivalent to not allowing the mayor to run for reelection)
                    IF (flag_no_run == 1) THEN
                       IF (DBLE(auditf)*dum_steal > 0.5d0) yy = -1000d0
!!$                       IF (DBLE(auditf)*dum_steal > 0.5d0 .AND. myrank == 0) THEN
!!$                          write(*,'(a,9f20.10,3i5)') 'emaxr nor', emaxr-run_shock_t, emaxnor, emaxr, run_shock_t, fr_stolen, ability, ln_relqcons, rundec_netW, msaving_bf, ub_sav, auditf, ipastf
!!$                          write(*,*) ''
!!$                       END IF
                    END IF
                    ! If we want to implement a policy in which the mayor looses the election with probability 1 if he steals more than x (equivalent to not allowing the mayor to run for reelection)
                    IF (flag_no_run_x == 1) THEN
                       IF (DBLE(auditf)*dum_steal_x > 0.5d0) yy = -1000d0
                    END IF
                    
                 END IF
                 
                 DEALLOCATE(betart,betanort)
                 
                 !We save the variables for the last period
80               simout(i,tt,1) = DBLE(i)
                 simout(i,tt,2) = DBLE(j)
                 simout(i,tt,3) = DBLE(tt)
                 simout(i,tt,4) = fr_stolen
                 !newqconsf is in levels
                 simout(i,tt,5) = qconsf
                 simout(i,tt,6) = mwagef
                 simout(i,tt,7) = DBLE(auditf)
                 simout(i,tt,8) = elap_vec(el_ap)
                 !These two variables are saved later in the file
!!$                 simout(i,tt,9) = fundsf
!!$                 simout(i,tt,10) = zzprf/EXP(zzpr_shock_t)
                 simout(i,tt,11) = ability
                 
                 !This is only useful is we write out data
                 IF (flag_simulate_data == 1) THEN
                    simoutt(i,j,tt,1) = DBLE(i) 
                    simoutt(i,j,tt,2) = DBLE(j) 
                    simoutt(i,j,tt,3) = DBLE(tt) 
                    simoutt(i,j,tt,4) = fr_stolen
                    !newqconsf is in levels
                    simoutt(i,j,tt,5) = qconsf
                    simoutt(i,j,tt,6) = mwagef
                    simoutt(i,j,tt,7) = DBLE(auditf)
                    simoutt(i,j,tt,8) = elap_vec(el_ap)
                    !These two variables are saved later in the file
!!$                    simoutt(i,j,tt,9) = fundsf
!!$                    simoutt(i,j,tt,10) = zzprf/EXP(zzpr_shock_t)
                    simoutt(i,j,tt,11) = ability
                 END IF
                 
                 thief=0
                 
                 ! fr_stolen is for the previous period
                 IF(fr_stolen > l_fr_stolen) thief =1

                 ! This corresponds to what we have in probmayor: Phi(Xb) =  N_cdf(Xb <= epsilon); hence, we need yy = Xb - epsilon >= 0 (all cases in which -epsilon is small enough to be smaller than Xb
                 IF (outelect == 1) THEN
                    
                    !The incumbent stays in power
                    nterms = nterms + 1
                    mayorage_elec = mayorage

!!$                    IF (myrank == 0) write(*,'(a,3i5,7f20.10,1f40.10)') 'ability in 2',i,j,tt, ability, sigma_a, ability_norm(i,j,tt), mu_a, d_run, emaxnor, emaxr, yy
                 
                 ELSE

                    !The challenger is elected
                    auditf = 0
                    ipastf = 1
                    nterms = 1

                    ! If we simulated the base case without selection on observables in the first and second terms, then flag_no_selec_obs = 1 and flag_no_selec_obs_from_1st = 1, and we assume that all mayors have the average variables observed in the data that are not mayor specific; if flag_no_selec_obs_from_1st == 1 then we set all observables equal to the average or median
!!$                    IF (flag_no_selec_obs == 1 .AND. flag_no_selec_mun_obs_from_1st == 1) THEN
!!$                       IF (flag_no_selec_obs_from_1st == 1) THEN


!!$                    write(*,'(a,6i4,2f20.10)') 'outelect_base 1', j, i, tt, flag_no_selec_obs, mayorage, ieducf,outelect_base, msavingf
                    
                    IF (imun == 1) THEN
                       
                       aenoc_mun = INT(uud(4) * (DBLE(countNnterms_mun1(nterms))-1.0d0)) + 1                       
                       mayorage = INT(nonpardistnoc_mun1(aenoc_mun,nterms,1))
                       ieducf = INT(nonpardistnoc_mun1(aenoc_mun,nterms,2))
                       msavingf = nonpardistnoc_mun1(aenoc_mun,nterms,5)
                       
                    ELSE IF (imun == 2) THEN
                       
                       aenoc_mun = INT(uud(4) * (DBLE(countNnterms_mun2(nterms))-1.0d0)) + 1
                       mayorage = INT(nonpardistnoc_mun2(aenoc_mun,nterms,1))
                       ieducf = INT(nonpardistnoc_mun2(aenoc_mun,nterms,2))
                       msavingf = nonpardistnoc_mun2(aenoc_mun,nterms,5)
                       
                    ELSE IF (imun == 3) THEN
                       
                       aenoc_mun = INT(uud(4) * (DBLE(countNnterms_mun3(nterms))-1.0d0)) + 1                       
                       mayorage = INT(nonpardistnoc_mun3(aenoc_mun,nterms,1))
                       ieducf = INT(nonpardistnoc_mun3(aenoc_mun,nterms,2))
                       msavingf = nonpardistnoc_mun3(aenoc_mun,nterms,5)
                       
                    END IF
                    
                    !For the second term we observe in the data, we use mayor's age and education from the data
                    IF (tt == 2) THEN

                       mayorage = INT(initdata(i,15))
                       ieducf = INT(initdata(i,16))
                       !This is only useful when we study selection to have the same distributio of observables as in tt = 1
!!$                       aenoc_mun = INT(uud(4) * (DBLE(Nactmun)-1.0d0)) + 1
!!$                       mayorage = INT(initdata(aenoc_mun,3))
!!$                       ieducf = INT(initdata(aenoc_mun,4))
!!$                       msavingf = initdata(aenoc_mun,18)/norm

                    END IF
                    
                    IF (flag_no_selec_obs == 1) THEN
                       IF (flag_no_selec_age == 1) mayorage = incumb_age
                       IF (flag_no_selec_educ == 1) ieducf = incumb_educ
                       IF (flag_no_selec_sav == 1) msavingf = msaving_base
                    END IF
                    
                    mayorage_elec = MIN(Retage-1,mayorage)
                    
!!$                    write(*,'(a,6i4,2f20.10)') 'outelect_base 2', j, i, tt, mayorage, mayorage_elec, ieducf, outelect_base, msavingf

                    !We draw ability from a log-normal distribution
                    ability = EXP(sigma_a*ability_norm(i,j,tt)+mu_a)
                    
                    !If we fix ability at some level
                    IF (flag_fixed_ability == 1) ability = fixed_ability
                    
!!$                    IF (myrank == 0) write(*,'(a,3i5,7f20.10,1f40.10)') 'ability cha 1',i,j,tt, ability, sigma_a, ability_norm(i,j,tt), mu_a, d_run, emaxnor, emaxr, yy
                    
                    !We draw the electoral appeal type
                    ! 1 is low electoral appeal, 2 is high electoral appeal, pappeal(1) is the probability of a low electoral appeal type
                    el_ap = 1
                    IF (uud(6) >= pappeal(1)) THEN
                       el_ap = 2
                    END IF
                    
                 END IF

!!$                 ! If we simulated the base case without selection on observables in the second term, then flag_no_selec_obs = 1 and we do not update the municipality's observables 
!!$                 IF (flag_no_selec_obs == 0) THEN
                 !For the second term we observe in the data, we use funds and the private inputs from the data
                 !In tt == 2 a mayor can never be in the third term because in tt == 1 they all start as first term mayors
                 IF (tt == 2) THEN
                    
                    fundsf = per_funds*initdata(i,14)/norm
                    zzprf = initdata(i,17)
                    
                 ELSE 
                    
                    IF (imun == 1) THEN
                       
                       aenoc_mun = INT(uud(4) * (DBLE(Ntotmunnoc_mun1)-1.0d0)) + 1
                       !Funds provided by central governement
                       fundsf = per_funds*nonpardistnoc_uncond_mun1(aenoc_mun,1)
                       
                       !Private gdp
                       zzprf = per_funds*nonpardistnoc_uncond_mun1(aenoc_mun,2)
                       
!!$                       write(*,'(a,1i6,3f18.8)') 'funds', Ntotmunnoc_mun1, uud(4), nonpardistnoc_uncond_mun1(aenoc_mun,1), nonpardistnoc_uncond_mun1(aenoc_mun,2)
                       
                    ELSE IF (imun == 2) THEN
                       
                       aenoc_mun = INT(uud(4) * (DBLE(Ntotmunnoc_mun2)-1.0d0)) + 1
                       !Funds provided by central governement
                       fundsf = per_funds*nonpardistnoc_uncond_mun2(aenoc_mun,1)
                       
                       !Private gdp
                       zzprf = per_funds*nonpardistnoc_uncond_mun2(aenoc_mun,2)
                       
!!$                       write(*,'(a,1i6,3f18.8)') 'funds', Ntotmunnoc_mun2, uud(4), nonpardistnoc_uncond_mun2(aenoc_mun,1), nonpardistnoc_uncond_mun2(aenoc_mun,2)
                       
                    ELSE IF (imun == 3) THEN
                       
                       aenoc_mun = INT(uud(4) * (DBLE(Ntotmunnoc_mun3)-1.0d0)) + 1
                       !Funds provided by central governement
                       fundsf = per_funds*nonpardistnoc_uncond_mun3(aenoc_mun,1)
                       
                       !Private gdp
                       zzprf = per_funds*nonpardistnoc_uncond_mun3(aenoc_mun,2)
                       
!!$                       write(*,'(a,1i6,3f18.8)') 'funds', Ntotmunnoc_mun3, uud(4), nonpardistnoc_uncond_mun3(aenoc_mun,1), nonpardistnoc_uncond_mun3(aenoc_mun,2)
                       
                    END IF
                 END IF
                
                 !WE SAVE funds AND private inputs HERE BECAUSE THE VARIABLES ARE DRAWN FROM A DISTRIBUTION THAT DEPENDS ON nterms WHICH IS DETERMINED AFTER THE ELECTIONS
                 simout(i,tt,9) = fundsf
!!$                 simout(i,tt,10) = zzprf/EXP(zzpr_shock_t)
                 simout(i,tt,10) = zzprf
                 
                 !This is only useful if we write out data
                 IF (flag_simulate_data == 1) THEN
                    simoutt(i,j,tt,9) = fundsf
!!$                    simoutt(i,j,tt,10) = zzprf/EXP(zzpr_shock_t)
                    simoutt(i,j,tt,10) = zzprf
                 END IF
                 
              END IF  !END IF: Time dependent stuff!

              IF (mayorage_elec < Retage -1 .AND. nterms < Nnterms) THEN
                 ALLOCATE(betat1(MMnew),betat2(MMnew),betat3(MMnew))
                 !inside mayorsolution.f90 we need the betas for t+1 for mayors with different types of stealing behavior the current period
                 IF (Npast > 1) THEN
                    betat1(:) = betanewcr(mayorage_elec+1,imun,el_ap,1,nterms+1,ipauf,ieducf,:)
                    betat2(:) = betanewcr(mayorage_elec+1,imun,el_ap,2,nterms+1,ipauf,ieducf,:)
                    betat3(:) = betanewcr(mayorage_elec+1,imun,el_ap,3,nterms+1,ipauf,ieducf,:)
                 ELSE
                    betat1(:) = betanewcr(mayorage_elec+1,imun,el_ap,1,nterms+1,ipauf,ieducf,:)
                 END IF
                 MMg = MMnew

!!$                 IF (el_ap == 1) THEN
!!$                    betat1(Nwealth+1:MMnew,:) = 0d0
!!$                    betat2(Nwealth+1:MMnew,:) = 0d0
!!$                    betat3(Nwealth+1:MMnew,:) = 0d0
!!$                 END IF

              ELSE IF(mayorage_elec == Retage -1 .OR. nterms >=Nnterms) THEN 
                 ALLOCATE(betat1(Nwealth),betat2(Nwealth),betat3(Nwealth))
                 !If they cannot run, it is indifferent whether they steal in the current period
                 IF (Npast > 1) THEN
                    betat1(:) = betanornew(mayorage_elec+1,imun,1,:,ipauf,ieducf)
                    betat2(:) = betanornew(mayorage_elec+1,imun,2,:,ipauf,ieducf)
                    betat3(:) = betanornew(mayorage_elec+1,imun,3,:,ipauf,ieducf)
                 ELSE
                    betat1(:) = betanornew(mayorage_elec+1,imun,1,:,ipauf,ieducf)
                 END IF
                 MMg = Nwealth
              END IF
!!$              IF (j == 21 .AND. i >= 3324 .AND. i < 3334 .AND. tt <= 3) flag_here = 178 
!!$              IF (j == 2 .AND. i == 3357 .AND. tt == 1 .AND. mayorage_elec == 5 .AND. imun == 3 .AND. ieducf == 2 .AND. ipauf == 1 .AND. ipastf == 1 .AND. is == 0 .AND. ifu == 0 .AND. ipr == 0 .AND. nterms == 1) flag_here = 456

              !We assign the probability of being audited, which depends on the audit regime ipau
              !It is used in mayorsolution0.f90 and mayorsolution.f90
              paudit_g = paudit(0,ipauf)
              IF (flag_audit_nterm == 1 .AND. (nterms == Nnterms .OR. mayorage_elec >= Retage-1)) paudit_g = paudit_nterm(0)

!!$              IF (i == 1 .AND. (j >= 1 .AND. j <= 5) .AND. tt == 1) THEN
!!$                 flag_here = 246
!!$                 ieducf = 2
!!$              END IF
!!$              IF(flag_here==246) write(*,'(a,12i5,6f18.9)') 'mayor info',i,j,tt,mayorage_elec,imun,ieducf,ipauf,ipastf,is,ifu,ipr,nterms,ability,mwagef,msavingf,fundsf,zzprf,elap_vec(el_ap)
!!$              IF(flag_here==246) write(*,'(a,3i4,15f50.45)') 'betas',i,j,tt,betat1(:,1)
!!$              IF(flag_here==246) write(*,'(a,3i4,15f50.45)') 'betas',i,j,tt,betat1(:,2)
!!$              IF(flag_here==246) write(*,'(a,3i4,15f50.45)') 'betas',i,j,tt,betat2(:,1)
!!$              IF(flag_here==246) write(*,'(a,3i4,15f50.45)') 'betas',i,j,tt,betat2(:,2)
!!$              IF(flag_here==246) write(*,'(a,3i4,15f50.45)') 'betas',i,j,tt,betat3(:,1)
!!$              IF(flag_here==246) write(*,'(a,3i4,15f50.45)') 'betas',i,j,tt,betat3(:,2)

!!$              write(*,'(a,6i4,2f20.10)') 'outelect_base 3', j, i, tt, mayorage, mayorage_elec, ieducf, outelect_base, msavingf

              ! We have to use the frac_fine drawn in this term (current)
              frac_fine_g = frac_fine_c
              CALL mayordecision(mayorage_elec,imun,ieducf,ability,ipauf,ipastf,mwagef,is,msavingf,ifu,fundsf,ipr,zzprf,nterms,newpconsf,newqconsf,zzpuf,mnewsavingf,newstealf,utilf,vemaxpf)

!!$              write(*,'(a,6i4,2f20.10)') 'outelect_base 4', j, i, tt, mayorage, mayorage_elec, ieducf, outelect_base, msavingf

!!$              flag_here = 0 

!!$              IF (myrank == 0) THEN
!!$                 IF (j == 50 .AND. i == 165 .AND. tt == 3) flag_here = 123
!!$                 !zzpuf is in per-capita terms
!!$                 !Ability type
!!$                 ability_g = ability
!!$                 CALL prodfunc(imun,ieducf,zzpuf,zzprf,qcons_t)
!!$                 flag_here = 0
!!$              END IF
!!$                 IF(flag_here==123) write(*,'(a,16f16.8)') 'qcons, prodfun, term 2',newqconsf,qcons_t*popsizeg,qcons_t,(zzpuf*norm)**alphaq(1)*zzprf**alphaq(2)*EXP(alphaq(3)*DBLE(ieducf)+alphaq(4))*ability,(zzpuf*norm)**alphaq(1)*zzprf**alphaq(2)*EXP(alphaq(3)*DBLE(ieducf)+alphaq(5))*ability,zzpuf,zzpuf*norm,zzprf,DBLE(ieducf),ability,alphaq(1:5),DBLE(nterms)
!!$                 IF(flag_here==123) WRITE(*,*) ''
!!$                 flag_here = 0
!!$              END IF

              ! We compute the value of the policy
              IF (flag_simulate_data == 1) THEN
                 ! We should compute the value function for the median resident only in the first period we consider
                 ! Since we compute the value function, it should include the effect of the policy for the rest of the median resident's life
                 ! THIS IS NOT TRUE BECAUSE IN tt == 1 ALL MAYORS ARE FIRST TERM AND THEY MAKE DIFFERNET DECISIONS
                 ! We use median age, education, and wealth conditional on the size of the municipality
                 med_edu_i = med_edu(imun)
                 med_ipast_i = med_ipast(imun)
                 med_age_g = med_age(imun)
                 med_wealth_g = med_wealth(imun)
                 pmwage_g = wage_pastmayor(med_age_g,imun,med_edu_i,med_ipast_i)
                 ELNqcons = newqconsf/((popsizeg)**eta)
                 ALLOCATE(betat(Nwealth))
                 MMg = Nwealth
                 betat(:) = betapmayornor(med_age_g+1,imun,med_edu_i,ipauf,med_ipast_i,:)
                 IF (flag_base == 1) THEN
                    
                    CALL notmayordecision(med_age_g,pmwage_g,med_wealth_g,ELNqcons)
                    ! We get the average value function for the municipality where the average is over simulations 
                    notmayor_value_t(i,j,tt) = vemaxoptpmg
                    ELNqcons_base_vec_t(i,j,tt) = ELNqcons
                       
!!$                    IF (myrank >= 0) write(*,'(a,8i5,5f18.8)') 'notmayor_value_t a', i,j,tt,imun,nterms,med_age_g,med_edu_i,med_ipast_i,pmwage_g,med_wealth_g,ELNqcons_g,notmayor_value_t(i,j,tt),vemaxoptpmg
!!$                    IF (myrank == 0) write(*,*) 'betat', betat(:)

                 ELSE

!!$                    IF (i >= 1 .AND. i <= 4 .AND. j == 1 .AND. tt >= 1) flag_fafb = 1

                    !We assign the citizens' value in the base case
                    notmayor_base_value = notmayor_value(i,j,tt)
                    ELNqcons_base = ELNqcons_base_vec(i,j,tt)

596                 CALL notmayordecision(med_age_g,pmwage_g,med_wealth_g,ELNqcons_base)
                    vemax_base = vemaxoptpmg

                    CALL notmayordecision(med_age_g,pmwage_g,med_wealth_g,ELNqcons)
                    vemax_policy = vemaxoptpmg

                    IF (vemax_base < vemax_policy) THEN
                       !Value function for the case with lower value, the value we want to achieve by lowering wealth
                       notmayor_low_value_g = vemax_base
                       !Value function for the other case
                       vemax_ini = vemax_policy
                       !Public consumption for the case with higher value, the public consumption we use when we compute the value function at lower values of wealth
                       ELNqcons_high_g = ELNqcons
                    ELSE
                       !Value function for the case with lower value, the value we want to achieve by lowering wealth
                       notmayor_low_value_g = vemax_policy
                       !Value function for the other case
                       vemax_ini = vemax_base
                       !Public consumption for the case with higher value, the public consumption we use when we compute the value function at lower values of wealth
                       ELNqcons_high_g = ELNqcons_base
                    END IF

!!$                       IF (flag_fafb == 1) write(*,'(A,2F20.10)') 'notmayor_base_value_g', notmayor_base_value_g, vemaxoptpmg
!!$                       IF (flag_fafb == 1) write(*,'(a,6i4,4f20.10)') 'notmayor_value_t a', i,j,med_age_g,imun,med_edu_i,med_ipast_i,pmwage_g,med_wealth_g,ELNqcons,ELNqcons_base

!!$                    write(*,'(A,2F20.10)') 'notmayor_base_value_g', notmayor_base_value_g, vemaxoptpmg
!!$                    IF (flag_fafb == 1) write(*,'(a,8i6)') 'notmayor_value_t dis', i,j,tt,imun,nterms,med_age_g,med_edu_i,med_ipast_i
!!$                    IF (flag_fafb == 1) write(*,'(a,7f30.20)') 'notmayor_value_t cont', pmwage_g,med_wealth_g,ELNqcons_base,ELNqcons,notmayor_base_value,vemax_base,vemax_ini
!!$                    IF (flag_fafb == 1) write(*,'(a,10f20.10)') 'betat', betat(:)
!!$                    IF (myrank >= 0) write(*,'(a,8i6)') 'notmayor_value_t dis', i,j,tt,imun,nterms,med_age_g,med_edu_i,med_ipast_i
!!$                    IF (myrank >= 0) write(*,'(a,7f30.20)') 'notmayor_value_t cont', pmwage_g,med_wealth_g,ELNqcons,ELNqcons_base,notmayor_base_low_g,vemax_base,vemax_policy
!!$                    IF (myrank >= 0) write(*,'(a,12i5,5f18.8)') 'mayor info',i,j,tt,mayorage_elec,imun,ieducf,ipauf,ipastf,is,ifu,ipr,nterms,ability,mwagef,msavingf,fundsf,zzprf

                    fr_wealth = 1d0
592                 a = med_wealth_g*fr_wealth
                    b = med_wealth_g
                    
                    fa = wealth_equivalent(a)
                    fb = wealth_equivalent(b)
                    
!!$                    IF (flag_fafb == 1) write(*,'(A,5F20.10)') 'fa, fb', a, fa, b, fb, fr_wealth
                    
                    IF(fa*fb <= 0.0d0) THEN
                       
                       !we look for the interior solution
                       wealth_star = zeroin(wealth_equivalent, a, b, 1d-10, 1d-10)
                       
                    ELSE
                       
                       IF (fb >= 0d0) THEN
                          fr_wealth = fr_wealth-0.2d0
                       ELSE
                          fr_wealth = fr_wealth+0.2d0
                       END IF
!!$                          flag_fafb = 1
!!$                          IF (myrank == 0) write(*,'(A,5F20.10)') 'WRONG STARTING POINTS', a, fa, b, fb, fr_wealth
!!$                          write(*,'(A,5F20.10)') 'WRONG STARTING POINTS', a, fa, b, fb, fr_wealth
                       GO TO 592
                       
                    END IF

                    IF (vemax_base < vemax_policy) THEN                    
                       will_to_pay = med_wealth_g - wealth_star
                    ELSE 
                       will_to_pay = wealth_star - med_wealth_g
                    END IF
                    av_will_to_pay_t(i,tt) = av_will_to_pay_t(i,tt) + will_to_pay

!!$                    IF (flag_fafb == 1) write(*,'(a,3i5,5f20.10)') 'will_to_pay', i,j,tt,will_to_pay, med_wealth_g, wealth_star, notmayor_low_value_g, vemax_ini
!!$                    IF (myrank >= 0) write(*,'(a,3i5,5f20.10)') 'will_to_pay', i,j,tt,will_to_pay, med_wealth_g, wealth_star, notmayor_base_value_g, vemax_ini
!!$                       IF (ABS(will_to_pay)>500d0) THEN
!!$                          flag_fafb = 1
!!$                          GO TO 596
!!$                       END IF
!!$
!!$                    flag_fafb = 0

                 END IF
                 DEALLOCATE(betat)
              END IF

              !We save optimal stealing without measurement errors
              newstealf_nome = newstealf
              !We add a measurement error to stealing: A municipality can make a clerical/administrative mistake that will be detected by the auditors
              IF (newstealf + mu_steal*fundsf + sigma_steal*fundsf*rnormprob(i,j,tt) > 0d0) THEN
                 newstealf = newstealf + mu_steal*fundsf + sigma_steal*fundsf*rnormprob(i,j,tt)
              ELSE
                 newstealf = 0d0
              END IF

!!$              IF (myrank == 0) write(*,'(a,6f16.10)') 'newsteal meas. errors after', newstealf, sigma_steal, fundsf, rnormprob(i,j,tt), sigma_steal*fundsf*rnormprob(i,j,tt), newstealf + sigma_steal*fundsf*rnormprob(i,j,tt)

              !If there is an infeasible solution we give the subsistence level to the mayor
              IF (utilf < check_bignumber .OR. ISNAN(mnewsavingf)) THEN 
                 infeasible = 1
                 WRITE(*,*) 'infeasible solution', newqconsf,mnewsavingf,newstealf,utilf,mayorage,imun,ieducf,ability,ist,auditf,mwagef,msavingf,fundsf,zzprf,yy,tt
                 newqconsf = subqcons
                 mnewsavingf = 0.0d0
                 newstealf = -999d0
                 utilf = subutil(INT(mayorage))
                 yy = -77777777777777d0
                 outelect = -777777
!!$                 flag_here = 890
!!$                 go to 890
              END IF

              IF (mayorage_elec < Retage) THEN
                 DEALLOCATE(betat1,betat2,betat3)
              END IF

              IF(tt == 3) THEN
                 !audit
!!$                 uud(5) = udraw(i,j,1,5)
                 pauditt(1:2) = paudit(0:1,ipau)
                 IF (flag_audit_nterm == 1 .AND. nterms == Nnterms) pauditt(1:2) = paudit_nterm(0:1)
                 CALL rmultinomial(uud(5), 2, pauditt, aud3)
                 auditf3 = aud3-1
!!$                 IF(newstealf > l_steal_prob/norm) thief3 = 1 
                 IF(newstealf/fundsf > l_fr_stolen) thief3 = 1 
              END IF

              !We save the variables related to the new mayor
              simout(i,tt,12) = msavingf
              simout(i,tt,13) = DBLE(ieducf)
              simout(i,tt,14) = ability
              simout(i,tt,15) = DBLE(outelect)
              simout(i,tt,16) = mnewsavingf 
              simout(i,tt,17) = newstealf
              !newqconsf is in levels
              simout(i,tt,18) = newqconsf
              simout(i,tt,19) = yy
              simout(i,tt,20) = DBLE(nterms)
              simout(i,tt,21) = zzpuf
              simout(i,tt,22) = DBLE(ipastf)
              simout(i,tt,23) = DBLE(imun)
              simout(i,tt,24) = will_to_pay
              simout(i,tt,25) = newstealf_nome
              simout(i,tt,26) = DBLE(auditf3)
              simout(i,tt,27) = DBLE(mayorage_elec)
              simout(i,tt,28) = popsize(imun)
              simout(i,tt,29) = stealf
              IF (tt > 1) THEN
                 simout(i,tt,30) = rnorm(i,j,tt)
              ELSE
                 simout(i,tt,30) = 0d0
              END IF
              simout(i,tt,31) = newpconsf
              simout(i,tt,32) = incumb_age
              simout(i,tt,33) = frac_fine_c
              simout(i,tt,34) = potential_wage
              simout(i,tt,35) = d_run
              simout(i,tt,36) = elap_vec(el_ap)

              !This is only useful if we write out data
              IF (flag_simulate_data == 1) THEN
                 simoutt(i,j,tt,12) = msavingf
                 simoutt(i,j,tt,13) = DBLE(ieducf)
                 simoutt(i,j,tt,14) = ability
                 simoutt(i,j,tt,15) = DBLE(outelect)
                 simoutt(i,j,tt,16) = mnewsavingf 
                 simoutt(i,j,tt,17) = newstealf
                 !newqconsf is in levels
                 simoutt(i,j,tt,18) = newqconsf
                 simoutt(i,j,tt,19) = yy
                 simoutt(i,j,tt,20) = DBLE(nterms)
                 simoutt(i,j,tt,21) = zzpuf
                 simoutt(i,j,tt,22) = DBLE(ipastf)
                 simoutt(i,j,tt,23) = DBLE(imun)
                 ! Extra for new moment: dummy for steal >0 for discrete level of stealing
                 simoutt(i,j,tt,24) = will_to_pay
                 simoutt(i,j,tt,25) = newstealf_nome
                 simoutt(i,j,tt,26) = DBLE(auditf3)
                 simoutt(i,j,tt,27) = DBLE(mayorage_elec)
                 simoutt(i,j,tt,28) = popsize(imun)
                 simoutt(i,j,tt,29) = stealf
                 IF (tt > 1) THEN
                    simoutt(i,j,tt,30) = rnorm(i,j,tt)
                 ELSE
                    simoutt(i,j,tt,30) = 0d0
                 END IF
                 simoutt(i,j,tt,31) = newpconsf
                 simoutt(i,j,tt,32) = incumb_age
                 simoutt(i,j,tt,33) = frac_fine_c
                 simoutt(i,j,tt,34) = potential_wage
                 simoutt(i,j,tt,35) = d_run
                 simoutt(i,j,tt,36) = elap_vec(el_ap)
              END IF

           END DO
           
        END DO

        ! We compute the numerator and denominator for all the moments

!!$        ! Moment for the identification of the relative taste for public consumption (we should use the ratio of average qcons to average savings
!!$        ! Average Public Consumption
!!$        numt(1) = numt(1) + SUM(simout(:,2:2,18)/simout(:,2:2,28), simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(simout(:,2:2,19) > -99999999999999d0))
!!$        dent(1) = dent(1) + 1d0
                
        ! Moments for the identification of fine parameter mu_a
        ! Average fraction of funds stolen conditional on second term mayor
        numt(1) = numt(1) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0))
        dent(1) = dent(1) + 1d0

        ! Average fraction of funds stolen conditional on second term mayor
        numt(2) = numt(2) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0))
        dent(2) = dent(2) + 1d0
                
        ! Moment to identify the fine parameter sigma_a (it is probably better to use the average of stealing squared).
        ! Now, Prob(Steal>0 & Auditf==1)/Prob(Audit==1)
        !Fraction of mayor who stole in the second term
!!$        numt(3) = numt(3) + DBLE(COUNT(simout(:,2:2,17)>l_fr_stolen .AND. INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0))/DBLE(COUNT(simout(:,2:2,17)>-1d0 .AND. INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0))
        !Fraction of mayor who stole
        numt(3) = numt(3) + DBLE(COUNT(simout(:,2:2,17)>l_fr_stolen .AND. simout(:,2:2,19) > -99999999999999d0))/DBLE(COUNT(simout(:,2:2,17)>-1d0 .AND. simout(:,2:2,19) > -99999999999999d0))
        dent(3) = dent(3) + 1d0

        ! Moment to identify the value of being in power:
        ! Fraction of mayors who did not steal and run for re-election; we should use the decision for period 3, which corresponds to the decision to run at the end of the term 2000-2004 in the data 
        thres_frac_stolen = l_fr_stolen
        ! We need to consider only mayors for which it was feasible to run (nterms < Nnterms, etc): simout(:,3:3,35) > -0.5d0
        numt(4) = numt(4) + DBLE(COUNT(simout(:,3:3,35) > 0.5d0 .AND. simout(:,3:3,7) > 0.5 .AND. simout(:,3:3,4) <= thres_frac_stolen .AND. simout(:,3:3,19) > -99999999999999d0))/DBLE(COUNT(simout(:,3:3,35) > -0.5d0 .AND. simout(:,3:3,7) > 0.5 .AND. simout(:,3:3,4) <= thres_frac_stolen .AND. simout(:,3:3,19) > -99999999999999d0))
!!$        numt(8) = numt(8) + DBLE(COUNT(simout(:,3:3,19) > -111111110d0 .AND. simout(:,3:3,4)*simout(:,3:3,7) <= thres_frac_stolen))/DBLE(COUNT(simout(:,3:3,19) > -111111112d0 .AND. simout(:,3:3,4)*simout(:,3:3,7) <= thres_frac_stolen))
        dent(4) = dent(4) + 1d0

!!$        if (myrank == 1) write(*,*) 'running', DBLE(COUNT(simout(:,3:3,35) > 0.5d0 .AND. simout(:,3:3,7) > 0.5 .AND. simout(:,3:3,4) <= thres_frac_stolen .AND. simout(:,3:3,19) > -99999999999999d0)), DBLE(COUNT(simout(:,3:3,35) > -0.5d0 .AND. simout(:,3:3,7) > 0.5 .AND. simout(:,3:3,4) <= thres_frac_stolen .AND. simout(:,3:3,19) > -99999999999999d0))

        !Probability of running for reelection conditional on stealing a positive amount of funds
        !We should use the decision for period 3, which corresponds to the decision to run at the end of the term 2000-2004 in the data
        ! We need to consider only mayors for which it was feasible to run (nterms < Nnterms, etc): simout(:,3:3,35) > -0.5d0
        numt(5) = numt(5) + DBLE(COUNT(simout(:,3:3,35) > 0.5d0 .AND. simout(:,3:3,7) > 0.5 .AND. simout(:,3:3,4) > thres_frac_stolen .AND. simout(:,3:3,19) > -99999999999999d0))/DBLE(COUNT(simout(:,3:3,35) > -0.5d0 .AND. simout(:,3:3,7) > 0.5 .AND. simout(:,3:3,4) > thres_frac_stolen .AND. simout(:,3:3,19) > -99999999999999d0))
!!$        numt(9) = numt(9) + DBLE(COUNT(simout(:,3:3,19) > -111111110d0 .AND. simout(:,3:3,4)*simout(:,3:3,7) > thres_frac_stolen))/DBLE(COUNT(simout(:,3:3,19) >= -111111112d0 .AND. simout(:,3:3,4)*simout(:,3:3,7) > thres_frac_stolen))
        dent(5) = dent(5) + 1d0

!!$        ! Moments to identify the constant in the production function and the mean ability
!!$        Nvarreg = 5
!!$        Ninst = (Nvarreg-1) + 1 ! Nvarreg-1: all regressors minus the endogenous variable; + 1: the instrument
!!$        ALLOCATE(y(Nactmun),x(Nactmun,Nvarreg),z(Nactmun,Ninst),betareg(Nvarreg),betaiv(Ninst),betatemp(Nactmun,Nvarreg),yhat(Nactmun))
!!$        y = 0d0
!!$        x = 0d0
!!$        z = 0d0
!!$        yhat = 0d0
!!$        betareg = 0d0
!!$        betaiv = 0d0
!!$        betatemp = 0d0
!!$
!!$        x(:,Nvarreg) = 1d0
!!$        ii = 0
!!$        tt = 2
!!$        DO imun = 1, Nactmun
!!$           ! We use both terms with a dummy for second term variables
!!$           IF (simout(imun,tt,18)/simout(imun,tt,28) > subqcons .AND. simout(imun,tt,21) > 0d0 .AND. simout(imun,tt,19) > -99999999999999d0) THEN
!!$
!!$              ii = ii + 1
!!$              ! qcons is in total terms, we need it in per-capita terms
!!$              y(ii) = LOG(simout(imun,tt,18)/simout(imun,tt,28))
!!$                            
!!$              ! zzpuf is already in per-capita terms
!!$              x(ii,1) = LOG(simout(imun,tt,21)*norm)
!!$              ! log zzprf
!!$              x(ii,2) = LOG(simout(imun,tt,10))
!!$              ! education
!!$              x(ii,3) = simout(imun,tt,13)
!!$              ! dummy for a second-term mayor
!!$              IF (simout(imun,tt,20) > 1.5d0) x(ii,4) = 1d0
!!$
!!$              z(ii,1:Nvarreg-1) = x(ii,2:Nvarreg)
!!$              !Funds
!!$              z(ii,Nvarreg-1+1) = LOG(simout(imun,tt,9)/simout(imun,tt,28))
!!$                                  
!!$           END IF
!!$        END DO
!!$
!!$        ! First Stage of the IV regression
!!$        CALL olsreg(ii,Ninst,x(1:ii,1),z(1:ii,:),betaiv,Vres,E2,out_err)
!!$
!!$        ! We replace the endogenous variable with its predicted values
!!$        x(1:ii,1) = 0d0
!!$        DO jj = 1, Ninst
!!$           x(1:ii,1) = x(1:ii,1) + betaiv(jj)*z(1:ii,jj)
!!$        END DO
!!$
!!$        ! Second Stage of the IV regression
!!$        ALLOCATE(temp_ols(ii,Nvarreg))
!!$        temp_ols(:,:) = x(1:ii,:)
!!$        CALL olsreg(ii,Nvarreg,y(1:ii),temp_ols(:,:),betareg,Vres,E2,out_err)
!!$        DEALLOCATE(temp_ols)
!!$
!!$         IF (out_err < 1d0) THEN
!!$
!!$            ! difference in constants estimated using d_2nd
!!$            numt(6) = numt(6) + betareg(Nvarreg-1)
!!$            dent(6) = dent(6) + 1d0
!!$            
!!$         END IF
!!$
!!$         DEALLOCATE(y,x,z,betareg,betaiv,yhat,betatemp)

!!$        newstealt = 0d0
!!$        !Median and 90th percentile
!!$        newstealt(:) = simout(:,2,17)/simout(:,2,9)
!!$
!!$        ! Median Stealing
!!$        ptest = 50d0
!!$        med_newsteal(1:1) = centile(newstealt(:), ptest)
!!$
!!$        numt(7) = numt(7) + med_newsteal(1)
!!$        dent(7) = dent(7) + 1d0

        !This moment identifies the prob of a high type
        !We compute the fraction of incumbents that are high appeal at the time of the election (simout(:,2:3,8) > 1.5d0) in this simulation
        !We should use the decision for period 3, which corresponds to the decision to run at the end of the term 2000-2004 in the data
        tt = 3
        numt(6) = numt(6) + DBLE(COUNT(simout(:,tt,8) > 0.5d0 .AND. simout(:,tt,35) > 0.5d0))/DBLE(COUNT(simout(:,tt,35) > 0.5d0))
        dent(6) = dent(6) + 1d0

        other_var_numt(1) = other_var_numt(1) + SUM(LOG(simout(:,2:2,14)), INT(simout(:,2:2,20))==1)/DBLE(COUNT(INT(simout(:,2:2,20))==1))
        other_var_dent(1) = other_var_dent(1) + 1d0

        other_var_numt(2) = other_var_numt(2) + SUM(LOG(simout(:,2:2,14)), INT(simout(:,2:2,20))==2)/DBLE(COUNT(INT(simout(:,2:2,20))==2))
        other_var_dent(2) = other_var_dent(2) + 1d0

        other_var_numt(3) = other_var_numt(3) + SUM(simout(:,2:2,14), INT(simout(:,2:2,20))==1)/DBLE(COUNT(INT(simout(:,2:2,20))==1))
        other_var_dent(3) = other_var_dent(3) + 1d0

        other_var_numt(4) = other_var_numt(4) + SUM(simout(:,2:2,14), INT(simout(:,2:2,20))==2)/DBLE(COUNT(INT(simout(:,2:2,20))==2))
        other_var_dent(4) = other_var_dent(4) + 1d0

        other_var_numt(5) = other_var_numt(5) + SUM(simout(:,2:2,18)/simout(:,2:2,28), simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(simout(:,2:2,19) > -99999999999999d0))
        other_var_dent(5) = other_var_dent(5) + 1d0

        other_var_numt(6) = other_var_numt(6) + SUM(simout(:,2:2,18)/simout(:,2:2,28), INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0))
        other_var_dent(6) = other_var_dent(6) + 1d0

        other_var_numt(7) = other_var_numt(7) + SUM(simout(:,2:2,18)/simout(:,2:2,28), INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0))
        other_var_dent(7) = other_var_dent(7) + 1d0

        other_var_numt(8) = other_var_numt(8) + SUM(simout(:,2:2,15), simout(:,2:2,35) > 0.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(simout(:,2:2,35) > 0.5d0 .AND. simout(:,2:2,19) > -99999999999999d0))
        other_var_dent(8) = other_var_dent(8) + 1d0

!!$
!!$        other_var_numt(1) = other_var_numt(1) + SUM(simout(:,2:3,18)/simout(:,2:3,28))/DBLE(2*Nactmun)
!!$        other_var_dent(1) = other_var_dent(1) + 1d0
!!$
!!$        other_var_numt(2) = other_var_numt(2) + SUM(simout(:,2:3,16))/DBLE(2*Nactmun)
!!$        other_var_dent(2) = other_var_dent(2) + 1d0

!!$        numt(10) = numt(10) + SUM(simout(:,2:2,14), INT(simout(:,2:2,20))==1)/DBLE(COUNT(INT(simout(:,2:2,20))==1))
!!$        dent(10) = dent(10) + 1d0
!!$
!!$        numt(11) = numt(11) + SUM(simout(:,2:2,14), INT(simout(:,2:2,20))==2)/DBLE(COUNT(INT(simout(:,2:2,20))==2))
!!$        dent(11) = dent(11) + 1d0

!!$        DO ii = 1, Nnterms            
!!$           DO jj = 1, Ntypes
!!$              frac_type_numt(jj,ii) = frac_type_numt(jj,ii)+DBLE(COUNT(INT(simout(:,2:2,14))==jj .AND. INT(simout(:,2:2,20))==ii))
!!$           END DO
!!$           frac_type_dent(ii) = frac_type_dent(ii)+DBLE(COUNT(INT(simout(:,2:2,20))==ii))
!!$        END DO
!!$
!!$        DO jj = 1, Ntypes
!!$           frac_type_run_numt(jj) = frac_type_run_numt(jj)+DBLE(COUNT(INT(simout(:,2:2,14))==jj .AND. simout(:,2:2,19) > -111111110d0))
!!$        END DO
!!$        frac_type_run_dent(1) = frac_type_run_dent(1)+DBLE(COUNT(simout(:,2:2,19) > -111111110d0))
!!$
!!$        DO ii = 1, Nnterms            
!!$           DO jj = 1, Nrcont
!!$              frac_el_app_numt(jj,ii) = frac_el_app_numt(jj,ii)+DBLE(COUNT(INT(simout(:,2:2,36))==INT(rcont(jj)) .AND. INT(simout(:,2:2,20))==ii))
!!$           END DO
!!$           frac_el_app_dent(ii) = frac_el_app_dent(ii)+DBLE(COUNT(INT(simout(:,2:2,20))==ii))
!!$        END DO
!!$
!!$        DO ii = 1, Nnterms            
!!$           DO jj = 1, Ntypes
!!$              DO el_ap = 1, Nelap
!!$                 frac_type_el_app_numt(jj,el_ap,ii) = frac_type_el_app_numt(jj,el_ap,ii)+DBLE(COUNT(INT(simout(:,2:2,14))==jj .AND. INT(simout(:,2:2,36))==INT(rcont(el_ap)) .AND. INT(simout(:,2:2,20))==ii))
!!$              END DO
!!$           END DO
!!$           frac_type_el_app_dent(ii) = frac_type_el_app_dent(ii)+DBLE(COUNT(INT(simout(:,2:2,20))==ii))
!!$        END DO
!!$
!!$        ! Prob(Steal>0 & Auditf==1)/Prob(Audit==1) in first term
!!$        other_var_numt(1) = other_var_numt(1) + (DBLE(COUNT(INT(simout(:,3,7))==1 .AND. INT(simout(:,3,24))==1 .AND. INT(simout(:,3,20))==1)) + DBLE(COUNT(INT(simout(:,3,26))==1 .AND. INT(simout(:,3,25))==1 .AND. INT(simout(:,3,20))==1)))/(DBLE(COUNT(INT(simout(:,3,7))==1 .AND. INT(simout(:,3,20))==1)) +  DBLE(COUNT(INT(simout(:,3,26))==1 .AND. INT(simout(:,3,20))==1)))
!!$        
!!$        other_var_dent(1) = other_var_dent(1) + 1d0
!!$        
!!$        ! Prob(Steal>0 & Auditf==1)/Prob(Audit==1) in second term
!!$        other_var_numt(2) = other_var_numt(2) + (DBLE(COUNT(INT(simout(:,3,7))==1 .AND. INT(simout(:,3,24))==1 .AND. INT(simout(:,3,20))==2)) + DBLE(COUNT(INT(simout(:,3,26))==1 .AND. INT(simout(:,3,25))==1 .AND. INT(simout(:,3,20))==2)))/(DBLE(COUNT(INT(simout(:,3,7))==1 .AND. INT(simout(:,3,20))==2)) +  DBLE(COUNT(INT(simout(:,3,26))==1 .AND. INT(simout(:,3,20))==2)))
!!$        
!!$        other_var_dent(2) = other_var_dent(2) + 1d0
!!$                
!!$        !Public good in second term if stealing > 0
!!$        other_var_numt(3) = other_var_numt(3) + SUM(simout(:,2:2,18), INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,17) > 0d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,17) > 0d0))
!!$        other_var_dent(3) = other_var_dent(3) + 1d0
!!$        
!!$        !Public good in second term if stealing == 0
!!$        other_var_numt(4) = other_var_numt(4) + SUM(simout(:,2:2,18), INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,17) <= 0d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,17) <= 0d0))
!!$        other_var_dent(4) = other_var_dent(4) + 1d0
!!$        
!!$        DO ii = 1, Ntypes
!!$           !Public good in first term if type = i
!!$           other_var_numt(4+ii) = other_var_numt(4+ii) + SUM(simout(:,2:2,18), INT(simout(:,2:2,20))==1 .AND. INT(simout(:,2:2,14))==ii)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. INT(simout(:,2:2,14))==ii))
!!$           other_var_dent(4+ii) = other_var_dent(4+ii) + 1d0
!!$           !Public good in second term if type = i
!!$           other_var_numt(4+Ntypes+ii) = other_var_numt(4+Ntypes+ii) + SUM(simout(:,2:2,18), INT(simout(:,2:2,20))==2 .AND. INT(simout(:,2:2,14))==ii)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. INT(simout(:,2:2,14))==ii))
!!$           other_var_dent(4+Ntypes+ii) = other_var_dent(4+Ntypes+ii) + 1d0
!!$           
!!$        END DO
!!$        
!!$        !Fraction of Mayors Not Running by Type
!!$        DO ii = 1, Ntypes
!!$           other_var_numt(4+2*Ntypes+ii) = other_var_numt(4+2*Ntypes+ii) + (DBLE(COUNT(INT(simout(:,2,19)) == -111111111 .AND. INT(simout(:,2,14))==ii)) + DBLE(COUNT(INT(simout(:,3,19)) == -111111111 .AND. INT(simout(:,3,14))==ii)))/(DBLE(COUNT(INT(simout(:,2,19)) >= -111111111 .AND. INT(simout(:,2,14))==ii)) + DBLE(COUNT(INT(simout(:,3,19)) >= -111111111 .AND. INT(simout(:,3,14))==ii)))
!!$           other_var_dent(4+2*Ntypes+ii) = other_var_dent(4+2*Ntypes+ii) + 1d0
!!$        END DO
!!$        
!!$        !Fraction of Mayors Not Running by Type if Audited
!!$        DO ii = 1, Ntypes
!!$           other_var_numt(4+3*Ntypes+ii) = other_var_numt(4+3*Ntypes+ii) + (DBLE(COUNT(INT(simout(:,2,19)) == -111111111 .AND. INT(simout(:,2,14))==ii .AND. INT(simout(:,2,7))==1)) + DBLE(COUNT(INT(simout(:,3,19)) == -111111111 .AND. INT(simout(:,3,14))==ii .AND. INT(simout(:,3,7))==1)))/(DBLE(COUNT(INT(simout(:,2,19)) >= -111111111 .AND. INT(simout(:,2,14))==ii .AND. INT(simout(:,2,7))==1)) + DBLE(COUNT(INT(simout(:,3,19)) >= -111111111 .AND. INT(simout(:,3,14))==ii .AND. INT(simout(:,3,7))==1)))
!!$           other_var_dent(4+3*Ntypes+ii) = other_var_dent(4+3*Ntypes+ii) + 1d0
!!$        END DO
!!$        
!!$        !Fraction of Mayors Not Running by Type if Audited and Caught Stealing
!!$        DO ii = 1, Ntypes
!!$           other_var_numt(4+4*Ntypes+ii) = other_var_numt(4+4*Ntypes+ii) + (DBLE(COUNT(INT(simout(:,2,19)) == -111111111 .AND. INT(simout(:,2,14))==ii .AND. INT(simout(:,2,7))==1 .AND. INT(simout(:,2,24))==1)) + DBLE(COUNT(INT(simout(:,3,19)) == -111111111 .AND. INT(simout(:,3,14))==ii .AND. INT(simout(:,3,7))==1 .AND. INT(simout(:,3,24))==1)))/(DBLE(COUNT(INT(simout(:,2,19)) >= -111111111 .AND. INT(simout(:,2,14))==ii .AND. INT(simout(:,2,7))==1 .AND. INT(simout(:,2,24))==1)) + DBLE(COUNT(INT(simout(:,3,19)) >= -111111111 .AND. INT(simout(:,3,14))==ii .AND. INT(simout(:,3,7))==1 .AND. INT(simout(:,3,24))==1)))
!!$           other_var_dent(4+4*Ntypes+ii) = other_var_dent(4+4*Ntypes+ii) + 1d0
!!$        END DO
!!$        
!!$        !Fraction of Mayors Not Running by Type if Audited and Caught Stealing More Than X
!!$        thres_steal = 10d0
!!$        DO ii = 1, Ntypes
!!$           IF ((DBLE(COUNT(INT(simout(:,2,19)) >= -111111111 .AND. INT(simout(:,2,14))==ii .AND. INT(simout(:,2,7))==1 .AND. simout(:,2,4) > thres_steal)) + DBLE(COUNT(INT(simout(:,3,19)) >= -111111111 .AND. INT(simout(:,3,14))==ii .AND. INT(simout(:,3,7))==1 .AND. simout(:,3,4) > thres_steal))) > 0d0) THEN
!!$              other_var_numt(4+5*Ntypes+ii) = other_var_numt(4+5*Ntypes+ii) + (DBLE(COUNT(INT(simout(:,2,19)) == -111111111 .AND. INT(simout(:,2,14))==ii .AND. INT(simout(:,2,7))==1 .AND. simout(:,2,4) > thres_steal)) + DBLE(COUNT(INT(simout(:,3,19)) == -111111111 .AND. INT(simout(:,3,14))==ii .AND. INT(simout(:,3,7))==1 .AND. simout(:,3,4) > thres_steal)))/(DBLE(COUNT(INT(simout(:,2,19)) >= -111111111 .AND. INT(simout(:,2,14))==ii .AND. INT(simout(:,2,7))==1 .AND. simout(:,2,4) > thres_steal)) + DBLE(COUNT(INT(simout(:,3,19)) >= -111111111 .AND. INT(simout(:,3,14))==ii .AND. INT(simout(:,3,7))==1 .AND. simout(:,3,4) > thres_steal)))
!!$           ELSE
!!$              other_var_numt(4+5*Ntypes+ii) = other_var_numt(4+5*Ntypes+ii) + 0d0
!!$           END IF
!!$           other_var_dent(4+5*Ntypes+ii) = other_var_dent(4+5*Ntypes+ii) + 1d0
!!$        END DO
!!$        
!!$        !Amount Stolen by Type First Term
!!$        DO ii = 1, Ntypes
!!$           other_var_numt(4+6*Ntypes+ii) = other_var_numt(4+6*Ntypes+ii) + SUM(simout(:,2:2,17), INT(simout(:,2:2,20))==1 .AND. INT(simout(:,2:2,14))==ii)/DBLE(COUNT(INT(simout(:,2:2,20))==1.AND. INT(simout(:,2:2,14))==ii))
!!$           other_var_dent(4+6*Ntypes+ii) = other_var_dent(4+6*Ntypes+ii) + 1d0
!!$        END DO
!!$        
!!$        !Amount Stolen by Type Second Term
!!$        DO ii = 1, Ntypes
!!$           other_var_numt(4+7*Ntypes+ii) = other_var_numt(4+7*Ntypes+ii) + SUM(simout(:,2:2,17), INT(simout(:,2:2,20))==2 .AND. INT(simout(:,2:2,14))==ii)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. INT(simout(:,2:2,14))==ii))
!!$           other_var_dent(4+7*Ntypes+ii) = other_var_dent(4+7*Ntypes+ii) + 1d0
!!$        END DO
!!$
!!$        ! Probability of winning the elections
!!$        other_var_numt(4+8*Ntypes+1) = other_var_numt(4+8*Ntypes+1) + SUM(simout(:,3:3,15), simout(:,3:3,19) > -11111111d0)/DBLE(COUNT(simout(:,3:3,19) > -11111111d0))
!!$        
!!$        other_var_dent(4+8*Ntypes+1) = other_var_dent(4+8*Ntypes+1) + 1d0
        
     END IF

  END DO

  CLOSE(18)

  IF (flag_simulate_data == 1) THEN
     
     ALLOCATE(simouttt(Nactmun,Nsim,Nsimterms,Nvar))
     simouttt = 0d0
     
     CALL MPI_ALLREDUCE(simoutt(1,1,1,1), simouttt(1,1,1,1), Nactmun*Nsim*Nsimterms*Nvar, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)

     IF (flag_base == 1) THEN

        CALL MPI_ALLREDUCE(notmayor_value_t(1,1,1), notmayor_value(1,1,1), Nactmun*Nsim*Nsimterms, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
        DEALLOCATE(notmayor_value_t)

        CALL MPI_ALLREDUCE(ELNqcons_base_vec_t(1,1,1), ELNqcons_base_vec(1,1,1), Nactmun*Nsim*Nsimterms, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
        DEALLOCATE(ELNqcons_base_vec_t)

     END IF

     av_will_to_pay = 0d0
     CALL MPI_ALLREDUCE(av_will_to_pay_t(1,1), av_will_to_pay(1,1), Nactmun*Nsimterms, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
     av_will_to_pay = av_will_to_pay/DBLE(Nsim)

!!$     DO i = 1, Nactmun
!!$        DO tt = 1, Nsimterms 
!!$           IF (myrank == 0) write(*,'(a,2i6,1f20.10)') 'av_will_to_pay', i, tt, av_will_to_pay(i,tt)
!!$        END DO
!!$     END DO
!!$
!!$     IF (myrank == 0) write(*,'(a,1f20.10)') 'av_will_to_pay over municipalities', SUM(av_will_to_pay(:,:))/DBLE(Nactmun*Nsimterms)

     IF (myrank == 0) THEN

        ! ESTIMATION OUTPUT FILE
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_estimation_Mar2021.dat", STATUS='REPLACE')
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_estimation_nf_Sep222017.dat", STATUS='REPLACE')
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_estimation_nse_Sep222017.dat", STATUS='REPLACE')
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_estimation_samefunds_Oct2017.dat", STATUS='REPLACE')

        ! BASE CASE
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_base_Mar2021.dat", STATUS='REPLACE')

        ! CONTROLLING FOR HETEROGENEITY
        !Output file when we simulate the policy that imposes a 1 term limit
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_1term_Mar2021.dat", STATUS='REPLACE') 
        !Output file when we simulate the policy that imposes a 1 term limit and we only remove selection on savings
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_1term_no_sel_sav_Mar2021.dat", STATUS='REPLACE')

        ! POLICIES OUTPUT FILE
        !INDIVIDUAL POLICIES
        !Output file when we simulate the policy that increases the audit probability to 0.168 (Lula's policy)
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_prau_Mar2021.dat", STATUS='REPLACE')
        !Output file when we simulate the policy that sets the probability of releelection to zero if someone is caught stealing (CRA policy)
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_norun_Mar2021.dat", STATUS='REPLACE')
        !Output file when we simulate a policy in which mayors who stole in the past are audited with probability 1
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_stealaudit_Mar2021.dat", STATUS='REPLACE')
        !Output file when we simulate the policy that increases the term limit to 3
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_terms_Mar2021.dat", STATUS='REPLACE')
        !Output file when we simulate the policy that doubles wages
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_wages_Mar2021.dat", STATUS='REPLACE')

        ! POLICY BUNDLES
        !Output file when we simulate the policy that increases the term limit to 3 and sets the probability of releelection to zero if someone is caught stealing (3 TERMS PLUS CRA)
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_terms_norun_Mar2021.dat", STATUS='REPLACE')
        !Output file when we simulate the three-term policy with higher probability of audit in the last term
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_terms_3termhigher_Mar2021.dat", STATUS='REPLACE')
        !Output file when we simulate the three-term policy plus no run with higher probability of audit in the last term
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_terms_norun_3termhigher_Mar2021.dat", STATUS='REPLACE')
        !Output file when we simulate the policy that increases the audit probability to 0.168 (Lula's policy) jointly with the audited-if-caught policy 
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_prau_stealaudit_Mar2021.dat", STATUS='REPLACE')
        !Output file when we simulate the no-run policy with higher probability of audit in the second term
!!$        OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/simdata_norun_2termhigher_Mar2021.dat", STATUS='REPLACE')

        DO j = 1, Nsim
           DO i = 1, Nactmun
              DO tt = 1, Nsimterms
                 WRITE(28,*) simouttt(i,j,tt,:)
              END DO
           END DO
        ENDDO
        CLOSE(28)

        IF (flag_base == 1) THEN      

           IF (myrank == 0) THEN
              OPEN(unit=28, FILE="~/project/MyPapers/Mayors/Estimation/Stata/Values_Base.dat", STATUS='REPLACE')
              DO j = 1, Nsim
                 DO i = 1, Nactmun
                    temp(1:Nsimterms) = notmayor_value(i,j,:)
                    temp(Nsimterms+1:2*Nsimterms) = ELNqcons_base_vec(i,j,:)
                    WRITE(28,*) temp
!!$                    IF (myrank == 0) write(*,'(a,2i6,6f20.10)') 'notmayor_value(i,j,:), ELNqcons_base_vec(i,j,:)', j, i, notmayor_value(i,j,:), ELNqcons_base_vec(i,j,:)
                 END DO
              END DO
           END IF
           CLOSE(28)
           
        END IF

     END IF

     CALL MPI_BARRIER(MPI_COMM_WORLD, ierr)
     
     DEALLOCATE(simouttt)
     DEALLOCATE(notmayor_value, ELNqcons_base_vec)     
     
  END IF
  
  IF(myrank == 0) write(*,*) ' ------------------------------------------------------------------'
  IF(myrank == 0) write(*,*) '  Compute the simulated moments for the estimation                 '
  IF(myrank == 0) write(*,*) ' ------------------------------------------------------------------'
  
  CALL MPI_BARRIER(MPI_COMM_WORLD, ierr)

  CALL MPI_ALLREDUCE(numt(1), num(1), Nmomtemp, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
  CALL MPI_ALLREDUCE(dent(1), den(1), Nmomtemp, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
!!$  CALL MPI_ALLREDUCE(frac_type_numt(1,1), frac_type_num(1,1), Nnterms*Ntypes, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
!!$  CALL MPI_ALLREDUCE(frac_type_dent(1), frac_type_den(1), Nnterms, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
!!$  CALL MPI_ALLREDUCE(frac_type_run_numt(1), frac_type_run_num(1), Ntypes, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
!!$  CALL MPI_ALLREDUCE(frac_type_run_dent(1), frac_type_run_den(1), 1, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
!!$  CALL MPI_ALLREDUCE(frac_el_app_numt(1,1), frac_el_app_num(1,1), Nnterms*2, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
!!$  CALL MPI_ALLREDUCE(frac_el_app_dent(1), frac_el_app_den(1), Nnterms, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
!!$  CALL MPI_ALLREDUCE(frac_type_el_app_numt(1,1,1), frac_type_el_app_num(1,1,1), Nnterms*2*Ntypes, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
!!$  CALL MPI_ALLREDUCE(frac_type_el_app_dent(1), frac_type_el_app_den(1), Nnterms, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
  CALL MPI_ALLREDUCE(other_var_numt(1), other_var_num(1), N_other_var, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
  CALL MPI_ALLREDUCE(other_var_dent(1), other_var_den(1), N_other_var, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)

  ! Print Moments !
  DO i = 1, nmom
     
     IF (ABS(den(i)) > 0.0d0) THEN
        simmom(i,1) = num(i) / den(i)
     ELSE
        simmom(i,1) = 0.0d0
     END IF

     IF (ISNAN(simmom(i,1))) THEN
        simmom(i,1) = 100d0
        IF (myrank == 0) write(*,'(a,1i4,3f20.10)') 'simmom NaN', i, simmom(i,1), num(i), den(i)
     END IF

  END DO
    
  DO i = 1, N_other_var
     other_var(i) = other_var_num(i) / other_var_den(i)
  END DO

  i = 0
  transform = 0d0  
  !Average Stealing for First Term Mayors
  i = i + 1
  IF (flag_sim == 0) THEN
     datamom(i,1) = data(i,1)*norm_factor(i)
     simmom(i,1) = simmom(i,1)*norm_factor(i)
  ELSE
     !We use the moments from the simulated data
     datamom(i,1) = simdatamom(i,1)*norm_factor(i)
  END IF
  transform(i,i)= 1.0d0*norm_factor(i)

  !Average Stealing for Second Term Mayors
  i = i + 1
  IF (flag_sim == 0) THEN
     datamom(i,1) = data(i,1)*norm_factor(i)
     simmom(i,1) = simmom(i,1)*norm_factor(i)
  ELSE
     !We use the moments from the simulated data
     datamom(i,1) = simdatamom(i,1)*norm_factor(i)
  END IF
  transform(i,i)= 1.0d0*norm_factor(i)
    
  !PR(steal & audit)/pr(audit)
  !Fraction mayors stealing in the second term
  i = i + 1
  IF (flag_sim == 0) THEN
     datamom(i,1) = data(i,1)*norm_factor(i)
     simmom(i,1) = simmom(i,1)*norm_factor(i)
  ELSE
     !We use the moments from the simulated data
     datamom(i,1) = simdatamom(i,1)*norm_factor(i)
  END IF
  transform(i,i)= 1.0d0*norm_factor(i)
    
  !Fraction of mayors running for reelection conditional on not stealing
  i = i + 1
  IF (flag_sim == 0) THEN
     datamom(i,1) = data(i,1)*norm_factor(i)
     simmom(i,1) = simmom(i,1)*norm_factor(i)
  ELSE
     !We use the moments from the simulated data
     datamom(i,1) = simdatamom(i,1)*norm_factor(i)
  END IF
  transform(i,i)= 1.0d0*norm_factor(i)
  
  !Probability of Running Conditional on Stealing
  i = i + 1
  IF (flag_sim == 0) THEN
     datamom(i,1) = data(i,1)*norm_factor(i)
     simmom(i,1) = simmom(i,1)*norm_factor(i)
  ELSE
     !We use the moments from the simulated data
     datamom(i,1) = simdatamom(i,1)*norm_factor(i)
  END IF
  transform(i,i)= 1.0d0*norm_factor(i)

!!$  ! Differences in constants between the first and second term
!!$  i = i + 1
!!$  IF (flag_sim == 0) THEN
!!$     datamom(i,1) = data(i,1)*norm_factor(i)
!!$     simmom(i,1) = simmom(i,1)*norm_factor(i)
!!$  ELSE
!!$     !We use the moments from the simulated data
!!$     datamom(i,1) = simdatamom(i,1)*norm_factor(i)
!!$  END IF
!!$  transform(i,i)= 1.0d0*norm_factor(i)

!!$  ! Median Fraction Stolen
!!$  i = i + 1
!!$  IF (flag_sim == 0) THEN
!!$     datamom(i,1) = data(i,1)*norm_factor(i)
!!$     simmom(i,1) = simmom(i,1)*norm_factor(i)
!!$  ELSE
!!$     !We use the moments from the simulated data
!!$     datamom(i,1) = simdatamom(i,1)*norm_factor(i)
!!$  END IF
!!$  transform(i,i)= 1.0d0*norm_factor(i)

  ! Probability of a high electoral appeal type
  i = i + 1
  IF (flag_sim == 0) THEN
     datamom(i,1) = data(i,1)*norm_factor(i)
     simmom(i,1) = simmom(i,1)*norm_factor(i)
  ELSE
     !We use the moments from the simulated data
     datamom(i,1) = simdatamom(i,1)*norm_factor(i)
  END IF
  transform(i,i)= 1.0d0*norm_factor(i)

  IF (myrank == 0) THEN
     i = 1
     write(*,*) ''
     write(*,*) 'Average Stealing in First Term, SIMULATION', simmom(i,1)
     write(*,*) 'Average Stealing in First Term, DATA', datamom(i,1)
     i = i + 1
     write(*,*) ''
     write(*,*) 'Average Stealing in Second Term, SIMULATION', simmom(i,1)
     write(*,*) 'Average Stealing in Second Term, DATA', datamom(i,1)
     i = i + 1
     write(*,*) ''
     write(*,*) 'PR(steal & audit)/pr(audit), SIMULATION', simmom(i,1)
     write(*,*) 'PR(steal & audit)/pr(audit), DATA', datamom(i,1)
     i = i + 1
     write(*,*) ''
     write(*,*) 'Fraction of Mayors Running Conditional on Not Stealing, SIMULATION', simmom(i,1)
     write(*,*) 'Fraction of Mayors Running Conditional on Not Stealing, DATA', datamom(i,1)
     i = i + 1
     write(*,*) ''
     write(*,*) 'Probability of running Conditional on Stealing, SIMULATION', simmom(i,1)
     write(*,*) 'Probability of running Conditional on Stealing,       DATA', datamom(i,1)
!!$     i = i + 1
!!$     write(*,*) ''
!!$     write(*,*) 'Differences in Production Function Constants, SIMULATION', simmom(i,1)
!!$     write(*,*) 'Differences in Production Function Constants, DATA', datamom(i,1)
!!$     i = i + 1
!!$     write(*,*) ''
!!$     write(*,*) 'Median Fraction Stolen SIMULATION', simmom(i,1)
!!$     write(*,*) 'Median Fraction Stolen, DATA', datamom(i,1)
     i = i + 1
     write(*,*) ''
     write(*,*) 'Prob High Appeal, SIMULATION', simmom(i,1)
     write(*,*) 'Prob High Appeal, DATA', datamom(i,1)
     write(*,*) ''
     write(*,*) 'Average Log Ability 1st Term, SIMULATION', other_var(1)
     write(*,*) ''
     write(*,*) 'Average Log Ability 2nd Term, SIMULATION', other_var(2)
     write(*,*) ''
     write(*,*) 'Difference in Log Average Ability between 2nd and 1st Term, SIMULATION', other_var(2)-other_var(1)
     write(*,*) ''
     write(*,*) 'Average Ability 1st Term, SIMULATION', other_var(3)
     write(*,*) ''
     write(*,*) 'Average Ability 2nd Term, SIMULATION', other_var(4)
     write(*,*) ''
     write(*,*) 'Difference in Average Ability between 2nd and 1st Term, SIMULATION', other_var(4)-other_var(3)
     write(*,*) ''
     write(*,*) 'Per-capita public consumption, SIMULATION', other_var(5)
     write(*,*) ''
     write(*,*) 'Per-capita public consumption 1st term, SIMULATION', other_var(6)     
     write(*,*) ''
     write(*,*) 'Per-capita public consumption 2nd term, SIMULATION', other_var(7)
     write(*,*) ''
     write(*,*) 'Fraction Elected, SIMULATION', other_var(8)
     write(*,*) ''
  END IF
  
  IF (flag_simulate_data == 1) THEN

     DEALLOCATE(simoutt)
     CALL MPI_FINALIZE(ierr)
     STOP
     
  END IF

  DEALLOCATE(simout)
  flag_simulation = 0
  
  CALL cpu_time(ts2)
  IF (myrank == 0) write(*,*) 'time simulation', ts2-ts1

END SUBROUTINE simulatemun
   
FUNCTION wealth_equivalent(wealthin)
  USE commonvars
  IMPLICIT NONE
  REAL(8), INTENT(IN) :: wealthin
  REAL(8)             :: wealth_equivalent

  EXTERNAL            :: notmayordecision
  
  CALL notmayordecision(med_age_g,pmwage_g,wealthin,ELNqcons_high_g)

  wealth_equivalent = vemaxoptpmg - notmayor_low_value_g

!!$  IF (myrank == 0) write(*,'(a,4f20.10)') 'wealth_equivalent', wealth_equivalent, wealthin, vemaxoptpmg, notmayor_low_value_g
!!$  IF (flag_fafb == 1) write(*,'(a,4f20.10)') 'wealth_equivalent', wealth_equivalent, wealthin, vemaxoptpmg, notmayor_low_value_g

END FUNCTION wealth_equivalent
